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SUMMARY 

This document discusses formulations and algorithms implemented in the 
MHOST finite element program developed under NASA Lewis Research 
Center Contract NAS3-23697. The code utilizes a novel concept of the mixed 
iterative solution technique for the efficient three dimensional computations 
of turbine engine hot section components. We first survey the state-of-the-art 
finite element technology at the time when the project was started (1983) and 
the major contributions made by this effort in the advancement of the 
methodology since then. Also the formulations and methodologies developed 
under this project are summarized in the first introductory chapter. The second 
chapter is devoted to the general framework of variational formulations and 
solution algorithms derived from the mixed three-field Hu-Washizu 
principle. This formulation enables us to use nodal interpolation for 
coordinates, displacements, strains and stresses. Algorithmic description of the 
mixed iterative method includes variations for the quasi static, transient 
dynamic and buckling analyses. The global-local analysis procedure referred to 
as the subelement refinement is developed in the framework of the mixed 
iterative solution, of which the detail is presented also in the second chapter. 
The third chapter discusses the numerically integrated isoparametric 
elements implemented in this framework. Methods to filter certain 
components of strain and project the element discontinuous quantities to the 
nodes are developed for a family of linear elements. In the fourth section, we 
describe integration algorithms for linear and nonlinear finite element 
equations included in the MHOST program. An appendix is added to outline 
the utilities handling the algebraic systems of linear equations such as the 
profile solver and subspace iterations for eigenvalue extraction. 
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1 . INTRODUCTION 


Introductory Remarks 

The main objective of the MHOST finite element program development is to construct 
an efficient analysis package for three dimensional models of gas turbine engine hot 
section components incorporating various constitutive models describing the material 
response under extreme thermal and mechanical loading conditions. Another 
requirement was to represent quantities involved in the analysis at nodes. These are 
geometry (coordinates), deformation (displacement and deformation gradient), strain, 
stress and material parameters. 

To fulfill the requirement, a new finite element formulation and solution 
algorithm is developed based on the mixed variational formulations. This concept 
enables us to introduce independent interpolation functions for all the physical 
quantities involved in the analysis, and to develop nodal interpolation algorithms. It 
was found from the literature survey that there had not been any formulation from 
which the efficient numerical procedure could be derived satisfying all the requirement 
stated in the first paragraph. 

The most significant accomplishment in this finite element method develop- 
ment is that the augmented Hu-Washizu principle was constructed from which both 
element formulations and solution algorithms are derived. In the formulation and 
computer program developments for the MHOST system, other major ingredients are: 

(i) The library of simple and efficient isoparametric elements ; 

(ii) A set of nonlinear constitutive equations for modelling high temperature in- 
elastic responses of hot section components; 

(iii) Efficient solution algorithms for quasi-static, dynamic and eigenvalue analysis ; 

(iv) Friendly user interface in defining the finite element models and loading con- 
ditions; also for graphic post-processing and other applications such as the per- 
turbation analysis for probabilistic structural analysis. 

The outline of this computational procedure differs radically from the con- 
ventional displacement method when applied to the linear elasticity. For nonlinear 
applications, the method works almost identical manner as the standard numerical 
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process. Indeed, one of the important findings in this project from the theoretical 
/numerical point of view is that the mixed iterative solution concept provides us with 
an rational framework to derive and analyze the solution strategies for nonlinear finite 
elements. 


User Input 

Material model 



Structural Model 

1 



Initial Solution 

► DISPLACEMENT PRECONDITIONING ◄ 

Numerical Filtering 

NODAL STRAIN PROJECTION 

Accurate Strain Field 

► NODAL STRESS RECOVERY 

Mixed Residual Vector 

EQUILIBRIUM ITERATIONS 

Improvement of Solution Quality 


Figure 1.1 The Mixed Iterative Process. A Schematic Flow Chart. 

Figure 1.1 illustrates the fundamental concept of mixed iterative solution 
strategy constructed directly from this mixed variational principle. The structural 
model is fed into the numerical computation in a usual way by constructing the 
displacement stiffness equations as the first step. Then the nodal strain is calculated by 
a mean square projection process, which eliminates the numerical instability in 
strain/stress approximations such as spurious oscillations. The inelastic material 
models are incorporated in the nodal stress integration in a modular fashion. The 
resulting nodally interpolated stress field is more accurate in comparison with the 
conventional displacement method. This improvement is brought back into the 
displacement finite element equations through the equilibrium iteration loop. 
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Additional computation required for the mixed iterative solution of linear 
problems is insignificant involving only a few re-solutions for the same stiffness 
matrix with modified residual vector. The nodal strain projection is carried out 
utilizing the diagonalized projection operator not requiring extra matrix 
manipulations. No additional computation is needed when the method is used for 
nonlinear computations. The modern iterative solution technology based on the quasi- 
Newton update is incorporated to further improve the performance of this solution 
method. 

The system of finite element equations solved in the MHOST program is 
symbolically shown in Figure 1.2 with the equivalent mixed stiffness equations derived 
by the direct elimination. The construction and inversion of this mixed finite element 
equation is prohibitively expensive in terms of both the computing time and the 
memory requirement. Hence the iterative solution algorithms become attractive 
alternatives to overcome this difficulty and take advantages of the mixed finite element 
formulations. 
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Figure 1.2 The Mixed Finite Element Equations. 


The linear convergence of the basic mixed iterative method is illustrated in 
Figure 1.3a. The difference of displacement and mixed stiffness matrices generates the 
driving residual vector for the iteration even for linear elastic problems. The quasi- 
Newton type update algorithms, with an example being illustrated in Figure 1.3b, 
improve the convergence rate of the iterative solution significantly. 

In the MHOST program, other than the constant metric iteration scheme 
(somewhat similar to the modified Newton iteration) a family of modem iterative 
solution algorithms are implemented in order to improve the efficiency including; the 
optional line search, the conjugate gradient, the secant implementation of Davidon 
rank-one quasi-Newton, BFGS rank-two quasi-Newton update algorithms. We observe 
that the displacement stiffness gives steeper gradient in the load-deflection relation 
than the exact one, whereas the mixed stiffness lies in the flexible side. Therefore the 
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iterative solution algorithm starting from the displacement method gradually 
recovering the mixed solution is almost always stable and convergent. 

To zoom-in to the local stress concentrations due to the embedded singularities 
such as holes and cracks, a scheme for global-local analysis referred to as the sub- 
element scheme is developed in the MHOST program. Elements requiring more 
detailed information in the global finite element computations are subdivided into a 
finer predefined mesh pattern, and the separate finite element computations are 
performed in these sub-element regions. The result is brought back to the global 
analysis by virtue of mixed iterative solution concept through the residual force vector 
generated by the sub-element stress field. 

Stiffness of the Displacement Preconditioner 




a) b) 

Figure 1.3 Behavior of the Mixed Iterative Solution. A Linear Problem. 

The numerically integrated isoparametric elements are used in this work. To be 
able to increase the computing speed, only the linear basis functions are used for the 
global finite element computations. The element library consists of: two node beam 
element (type 98), four node plane stress element (type 3); plane strain element (type 
11); axisymmetric solid element (type 10); three dimensional solid element (type 7); 
and, shell element (type 75). To be able to improve the accuracy, the selectively reduced 
integration scheme is incorporated with necessary refinements. Those are a rational 
coordinate transformation and strain filtering mechanism based on the polar decom- 
position of the isoparametric transformation and a projection method (often referred 
to as the variational recovery method) from the element discontinuous strain appro- 
ximation to the nodal interpolation. A set of elements with high coarse mesh accuracy 
derived from the assumed stress concept (often referred to as the assumed stress hybrid 
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elements) is derived and included in the element library. Those are: plane stress (type 
151), plane strain (type 152), axisymmetric solid (type 153), and three-dimensional solid 
(type 154). 

The material models included in the MHOST program are: 

(i) The linear elasticity with optional temperature dependency and anisotropy ; 

(ii) The simplified plasticity (the secant elastic modulus is calculated in the stress 
recovery which is then used to generate the secant stiffness matrix representing 
the plastic response under the monotonic loading situations ); 

(iii) The von Mises plasticity with the temperature dependent yield stress and the 
anisotropic yield surface definitions being the optional features. The radial 
return algorithm is used for the integration of the elastic-plastic material 
response over the given strain history. The general background of this type of 
algorithms is documented by HUGHES (1984). The finite deformation version 
of this algorithm implemented in this code is documented in NAGTEGAAL 
(1985). The effects of creep are taken into account in an explicit time integration 
algorithm constructed inside of the quasi-static solution subsystem. 

(iv) The unified viscoplastic constitutive equation originally developed by 
WALKER (1982) and implemented in a general purpose finite element package 
by CASSENTI (1983) is implemented in MHOST as a part of the material library. 
A straightforward explicit time integration with the substepping is used as the 
driving mechanism for the stress recovery. The temperature dependent 
isotropic elastic modulus is used as the preconditioner of the iterative solution 
incorporating this material model. 


The mechanism to incorporate other material models are provided in this computer 
program. As shown in Figure 1.1, the material definition is treated as another interface 
of the mixed iterative finite element method to the real world. The concept of user 
subroutines makes a wide range of different constitutive models easily implemented. 


Historical Background 

The research directions in the improvement of performance and applicability of the 
finite element method have continued to be: 

(i) Mixed and hybrid formulations to improve the accuracy of the finite elements ; 

(ii) Global solution procedures for linear and nonlinear finite element equations to 
reduce both the memory requirement and number of arithmetic operations 
with particular emphasis on the iterative solution strategies ; and, 
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(iii) Post-processing methods not only to accurately estimate errors a posteriori 
leading to adaptive mesh refinement strategies but also to extract more 
information out of a finite element solution. 

Equipped with the modern numerical technology and the advanced computing 
machinery, the application of finite elements expands its limit. The use of large 
meshes representing fine detail of the engineering problems has become useful. 
Sophisticated numerical methods have been developed to deal with combined 
nonlinear kinematics and material responses in a highly efficient manner. 

Reports, books and journal articles related to the non- standard finite element 
methods are surveyed with the main objective being the search of useful numerical 
technology in the framework of mixed iterative solution algorithms. In the past, the 
non-standard finite element methods such as the mixed and hybrid method have not 
been exploited in a systematic manner. Few exceptions are the hybrid elements for 
plates and shell analysis which eventually generates the element stiffness equations 
and the mixed elements for incompressible problems based on the Lagrangian 
functional used by HERRMANN (1965). 

As documented in TODD, CASSENTI, NAKAZAWA, BANERJEE (1984) and the 
research papers by NAKAZAWA (1984A,B), ZIENKIEWICZ, VTLOTTE, NAKAZAWA, 
TOYOSHIMA (1984), ZIENKIEWICZ, LI, NAKAZAWA (1985), the iterative algorithms 
to develop continuous stress field such as developed by CANTIN, LOUBIGNAC, 
TOUZOT (1978) can be identified as the mixed finite element method. It is important to 
note that, in constructing the iterative methods, the use of Hu-Washizu variational 
principle is crucial to set up the practically useful algorithms, although a contrary 
remark is made in the classical textbook by ZIENKIEWICZ (1977) (page 310). In the 
recent paper by SIMO, TAYLOR, PISTER (1984), the same observation is made and a 
finite element method is constructed from the Hu-Washizu principle with application 
to the finite deformation plasticity. 

Except for a few primitive papers such as mentioned above and the literature 
which appeared after the present development effort was initiated, no directly relevant 
work has been published on construction and iterative solution for the mixed finite 
element equation derived from the Hu-Washizu form. There are, however, a number 
of papers indirectly relevant to and somehow useful for the further development of 
solution strategy employed in the MHOST code. 

The algorithm is viewed as composed of three steps. The linearized 
momentum equation is solved in terms of displacement vector for the pre- 
conditioning purposes. Then , the post- processing algorithm is entered to gene- 
rate the strain field based on the mixed interpolation which is the used for the 
integration of the constitutive equation. The third step is the equilibrium 
iteration to satisfy the nonlinear virtual work equation with respect to the stress 
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field interpolated in a mixed manner. 

The quality of displacement field generated by the preconditioning contributes 
significantly to the overall quality of the solution as well as the convergence chara- 
cteristics. The use of equivalent stiffness to the element discontinuous strain mixed 
forms is a possible numerical strategy for these purposes [MALKUS, HUGHES (1978), 
HUGHES, MALKUS (1983)]. In particular, application to plates and shells in recent 
developments for the lower order element technology indicate possible improvements 
for the preconditioning operations [HUGHES, TEZDUYAR (1981)] with further effi- 
ciency gained by using the lower order quadrature [BELYTSCHKO, TSAY (1983)]. 

The strain recovery algorithm based on the mixed interpolation is found to be 
virtually identical to the classical methodology based on the consistent conjugate stress 
distribution studied by ODEN, BRAUCHU (1971), ODEN, REDDY (1973). In the recent 
papers, BABUSKA, MILLER (1984) present a systematic method to construct and 
analyze the post-processing algorithms. It is claimed that there are postprocessing 
procedures for various quantities which have the same order of convergence rate as the 
energy error in the displacement finite element approximation. This statement agrees 
quite well with the experimental observation by OWA (1982), NAKAZAWA, OWA, 
ZIENKIEWICZ (1985). Also, the super convergent results in terms of stress and strain 
as well as the displacement observed in the first year exercise of this project may be due 
to the higher order convergence rate of the post-processing algorithm used for the 
strain recovery computation. No literature is available on the effect of numerical 
quadrature in the postprocessing algorithms except a recent technical note by SIMO, 
HUGHES (1985) on the assumed strain, so called B-bar type algorithms for continuum 
elements and plates and shells. 

The integration algorithms of the nonlinear constitutive equations, in particular 
rate independent plasticity, are well- established as indicated in the Task I literature 
survey [FAHYRE, HUGHES (1983)] and one of the most reliable algorithm based on the 
radial return concept [NAGTEGAAL, DE JONG (1981)] is implemented. No systematic 
effort has yet been reported of the algorithm development for the advanced viscoplastic 
constitutive models including the temperature effect. Further literature search and 
original investigation is required in this field. 

Methods of equilibrium iteration have been investigated in recent years and a 
number of useful papers and reports have appeared, some of which are mentioned in 
the previous literature survey report. A survey and series of experiment reported by 
ABDEL RAHMAN (1982) is found a useful collection of algorithms and numerical 
results with application to the nonlinear plates. Basic algorithms of Newton-Raphson 
and modified- and quasi-Newton methods are compared in the displacement method 
framework. As demonstrated later in this report, algorithms discussed in a classical 
report by MATTHIES, STRANG (1979) are found usable even in the content of the 
mixed iterative method. 

The algorithms and convergence arguments directly applicable to the present 
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framework are only available from the literature on the augmented Lagrangian 
methods for the quadratic minimization problem with linear equality constraints such 
as incompressibility and the Dirichlet boundary condition. Possible improvement of 
the iterative procedure is indicated in FORTIN, GLOWINSKI (1982) and the numerical 
test examples studied by NAKAZAWA, VILOTTE, ZIENKIEWICZ (1984), 
ZIENKIEWICZ, VILOTTE, NAKAZAWA, TOYOSHIMA (1984) shows the significant 
improvement obtained by using such algorithms for the analysis of Stokes' flow. The 
mathematical discussions and algorithms are also found in GIRAULT, RAVIART 
(1979), TEMAM (1977). The original idea of the mixed iterative solution is, indeed, 
found in a historical work by ARROW, HURWICZ, UZAWA (1968) and the class of 
iterative algorithms for the mixed problems is referred to as the Uzawa method. 

In the context of linear elastic finite element analysis, the use of mixed 
approximations and equilibrium iteration has appeared in an ad hoc fashion repeatedly 
in the history other than the work by CAN TIN, et al. (1979). For nonconforming plate 
bending elements, CRISFIELD (1975) has proposed an iterative algorithm to improve 
the property of the numerical process. In the paper, the similarity of the concept based 
on a modified Hellinger-Reissner principle to the initial strain method iteration is 
pointed out. Also, an important observation is that all the significant changes occur in 
the first iteration of the mixed process of plates. 

The implementation of mixed finite elements without iterative solution is the 
topic extensively studied several decades ago mainly with application to the linear 
problems in mechanics as discussed in a survey by ZIENKIEWICZ (1977). These 
experiences are found useful to avoid possible numerical difficulties particular in this 
methodology. This is an important issue to be investigated. The characteristics of the 
particular algorithm used as the post-processor needs to be understood in the mixed 
method framework so as to fulfill the necessary conditions for the stability and 
convergence. Regardless of the solution algorithm, the mixed method needs to satisfy 
the stability condition referred to as the Babuska-Brezzi condition in some sense 
[BABUSKA (1973), BREZZI (1974)). However, when the equal order interpolations of 
displacement and stress are used, possible dissatisfaction of this condition is indicated 
by ODEN (1982). Possible oscillations of the numerical solution are indicated therein 
under special circumstances. 

As experienced in the mixed /penalty finite element computations for 
incompressible problems, however, the implication of Babuska-Brezzi condition is not 
quite clear. As discussed by ZIENKIEWICZ, NAKAZAWA (1982), the stability can be 
achieved by using a class of unstable elements which violates the condition a priori but 
produces stable results incorporating a post-processing algorithm which satisfies the 
necessary condition for the stability 

Modem development of the mixed finite element methods is mostly to set up 
the displacement stiffness matrix from the element discontinuous approximation for 
additional variables [TAYLOR, ZIENKIEWICZ (1982)]. The derivation of this class of 
methods is based on the equivalence theorem stated in MALKUS, HUGHES (1978). An 
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important development along this line is a generalization of the equivalence theorem 
proposed in ZIENKIEWICZ, NAKAZAWA (1984) indicating that for all the nume- 
rically integrated displacement finite element method, there exists an equivalent class 
of mixed methods based on the Hu-Washizu principle. This observation provides 
some insights to the iterative process with particular application to inelastic problems. 
As mentioned earlier, these developments of reducible mixed forms at the element 
level are useful mainly to construct the preconditioning system of equations as well as 
the fine tuning of the equilibrium iteration. 

In developing high performance finite elements, the Hu-Washizu principle has 
become the most popular variational formulation from which numerical methods are 
derived. It seems the most natural approach, for instance, to derive algorithms for 
probabilistic structural analysis methods. Also, methods to evaluate elements based on 
mixed variational principles have been developed. In modern literature, the word 
mixed and hybrid finite elements are often used interchangeably and ambiguously. 
Seemingly, the mixed methods derived from piecewise discontinuous stress 
approximation are usually referred to as the hybrid formulations and those derived 
from independent strain interpolation and other mixed formulations are referred to as 
the mixed forms. It is interesting to note that the word hybrid is used by a certain 
sector of computational mechanics community referring to a class of methods com- 
bining the finite element and boundary integral formulations. 

Attention has recently been paid to the equal order interpolation for the mixed 
finite element processes. Series of papers and reports have been coming out from not 
only this project but also the research group at Swansea. A number of academic 
research projects have been initiated in the recent few years. 

The interest in the Hu-Washizu principle is rejuvenated by the untimely death 
of one of the original investigators and the publication of his memorial volume [PLAN 
(1984)]. Utilization of this principle for the analysis of existing methods such as 
numerically integrated displacement formulation (ZIENKIEWICZ, NAKAZAWA 
(1984)] appears in literature frequently. Belytschko (1986) reviews the technology 
available for the thick plate and shell element formulation based on the Reissner- 
Mindlin theory by this principle. 

Efforts to use the mixed Hu-Washizu finite element equations as a driving 
mechanism for the iterative solution of linear and nonlinear problems are reported in 
recent publications [NAKAZAWA, NAGTEGAAL (1986), NAKAZAWA, SPIEGEL 
(1986)]. With particular application to transient dynamics, a version of the mixed 
iterative solution algorithm combined with a second order, single step time integration 
operator is proposed and applied to a number of linear elastic problems by 
ZIENKIEWICZ, LI, NAKAZAWA (1986). 

An analysis of a certain class of mixed finite elements is discussed by LE TALLAC 
(1986). An iterative algorithm based on the Hu-Washizu principle is derived for rate- 
independent deviatoric plasticity. The elastic strain is used as a variable instead of the 
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total strain. This formal alteration of the formulation in conjunction with the radial 
return algorithm being hard-wired to the analysis enabled him to discuss the existence 
of solution which leads us to the possible formalization of convergence characteristics. 

For the evaluation of element formulations and the first step of finite element 
code validation, the patch test invented by Irons and documented in BAZELEY, 
CHEUNG, IRONS, ZIENKIEWICZ (1956) has been used by many researchers and 
commercial code developers. NAGTEGAAL, NAKAZAWA, TATEISHI (1986) discuss 
the development of an eight-node shell element based on the assumed strain approach 
and demonstrate its validity by performing an extensive set of patch test. 

The concept of patch test is re-examined in a recent paper by Taylor, Simo, 
ZIENKIEWICZ, CHAN (1986) where the equivalence between the satisfaction of the 
patch test and the consistency and stability of finite element approximations is 
discussed. Note that the argument by STUMMEL (1980) seriously questioning the 
validity of the patch test is countered in this paper. 

An important extension of the patch test is proposed in a more recent paper by 
ZIENKIEWICZ, QU, TAYLOR, NAKAZAWA (1986). Assuming that the consistency of 
the approximations is satisfied, the test focuses upon the stability of mixed finite 
elements implied by the Babuska-Brezzi condition [BABUSKA (1973), BREZZI (1974)]. 
The procedure is to count the minimum possible number of active degrees-of- freedom 
in a patch and compare it with the maximum possible number of constraint degrees-of- 
freedom. If no constraints are imposed on nodal stresses and strains, the mixed, 
iterative solution procedure with continuous interpolation passes the test and is 
assumed stable. An extension of this stability patch test to three field formulations of 
Hu-Washizu type and for plates and shells derived from the Reissner-Mindlin theory 
is being investigated by ZIENKIEWICZ (1986) and his co-workers. 

The relation between the stability patch test and the mathematically derived 
Babuska-Brezzi condition is also being actively investigated. 

For the numerical modeling of singularities embedded in the structures, papers 
on the fracture mechanics are surveyed. For the displacement finite element 
technology, the report by FINE (1984) covers the technology to date. Except for a few 
academic exercises such as reported in ANNIGERI (1984), not many papers and research 
reports are available to mixed and hybrid finite elements with application to fracture 
mechanics, in particular, nonlinear material problems. A paper by BABUSKA, MILLER 
(1984) on the post-processing to detect the stress intensity factor is found useful to 
construct and validate the numerical algorithm to be implemented in the MHOST 
program. The general strategy for the numerical post-processing is extended to deal 
with problems with singularities. 

A series of papers by BABUSKA et. al. (1984) discusses the possible utilization of 
the post-processing technology and the adaptive mesh refinement. The concept is 
closely related to the mixed iterative solution algorithms in conjunction with the 
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subelement calculation as discussed later in this report. Literature on the adaptivity and 
a posteriori (See GAG O (1982) and literature cited therein) error estimate is found 
useful in this line of development. 

The methodology to compute the approximation error in a finite element 
solution has been improved significantly in the last few years. Also, the adaptive mesh 
refinement algorithms based on the spatial distribution of error indicator are developed 
and applied to a wide range of linear and nonlinear problems. A book on this subject is 
compiled by Gago which represents the state-of-the-art [BABUSKA, ZIENKIEWICZ, 
GAGO, OLIVERLA (1986)]. 

The papers included in the book indicate that a posterior error estimate and 
adaptive refinement methodologies are well understood mathematically and tested for 
a wide range of linear elastic problems. However, the computational processes look 
overly complicated and extensions to nonlinear analysis of solids and structures will 
need further research and development. Application of the adaptive mesh refinement 
concept has been demonstrated tractable using a relatively simple error indicator with 
application to aerospace flow problems given by the compressible Euler equations. 
Note that the adaptive solution of compressible Euler equations are the first systematic 
attempt of extending the technology to nonlinear problems. The results included in 
the book shows the potential advantage of the concept over the commonly used finite 
difference technology. 

A simple algorithm to estimate error a posteriori is proposed by ZIENKIEWICZ, 
ZHU (1986). The logic is based on the difference of quantities obtained directly at the 
elements and the same information projected by nodal interpolation functions. The 
nodal projection techniques used in the mixed, iterative process to recover the nodal 
strain is the simple and effective way to calculate errors. This process is relatively 
simple to implement in comparison with seemingly complicated mathematically 
derived a posterior error estimate algorithms. Examples shown in the paper indicate 
the efficiency of the proposed process. However, rigorous mathematical analysis is still 
to be done to verify the methodology. It is interesting to note that the residual vector 
calculated in MHOST at the zeroth iteration of the equal order mixed iterative solution 
is indeed the Zienkiewicz-Zhu error indicator weight-averaged at nodes. 

Using triangular elements, efforts to combine automatic mesh-generation and 
adaptive refinement concepts have been reported. Examples included in 
ZIENKIEWICZ, ZHU (1986) uses this idea. In application of adaptive mesh refinement 
to compressible flow problems, a number of papers are published in which triangular 
mesh generation methods were discussed with local refinement strategy to capture 
shocks. PERAIRE, MORGAN, ZIENKIEWICZ (1986), LOHNER (1986) discuss adaptive 
mesh refinement and de-refinement algorithm designed for the triangles and speculate 
the applicability of the concept to three-dimensional computations. The use of 
automatic mesh generation of triangular and tetrahedral elements has been found a 
useful trick in finite difference flow computations using unstructured grid. JAMESON 
(1986) demonstrated possible utilization of this technology for geometrically complex 
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problems. 

In the area of a posteriori error estimates, post-processing and adaptive mesh 
refinement, useful technology has been developed for linear structural analysis and 
certain nonlinear problems. However, further research and development including 
extensive numerical experiments are required to establish a methodology usable for 
three- dimensional inelastic analyses involving complex loading histories. 

A series of note by AXELSSON (1976, 1978) indicates the possible utilization of 
successive relaxation techniques for the solution of finite element equations, in 
particular, the mixed system of equations. However, no evidence that such algorithms 
can be used for nonlinear solution is provided in those publications. 

Application of finite element methods to nonlinear structural analysis involves 
increasingly complicated geometrical and material models and utilizes performance of 
currently available supercomputing facilities to their limits. Codes specifically designed 
for supercomputers have demonstrated their potential. Papers have appeared recently 
discussing performance improvement of finite element codes on vector and parallel 
machines utilizing algorithms tailored for specific computer environment. BENSON, 
HALLQUIST (1986) discusses the rigid body algorithm implemented in DYNA code 
which reduced the amount of computations significantly in explicit finite element time 
integrations. Implementation of an iterative solution algorithm based on the element- 
by-element preconditioner in NIKE code is discussed by HUGHES, FERENCZ, .. 
HALLQUIST (1986). The utilization of supercomputers in an engineering 
environment involving the development of these codes is presented by GOUDREAU, 
BENSON, HALLQUIST, KAY, ROSINSKY, SACKETT (1986). 

Investigations of vectorizability of basic finite element operations continue. For 
instance, the detail of vector-matrix multiply in element-by-element manner is 
discussed including timing figures for various finite element models by HAYES, 
DEVLOO (1986). An algorithm to potentially take full advantage of vector machines to 
solve a system of ordinary differential equations is proposed by Brown, Hindmarsh 
(1986). Iterative algorithms based on the ORTHOMIN method and incomplete 
factorization are studied by ZYVOLOSK3 (1986). It is anticipated that the development 
of algorithms closely tied to the actual data manipulations and arithmetic operations in 
computing machineries will eventually lead us to write faster finite element codes not 
only on supercomputers but also on smaller scalar machines. 

Utilization of parallel computing for the equation solution is a subject 
investigated extensivelyTn the last few years. A series of papers by Utku, Melosh and 
their collaborators look into the possible parallelization of direct stiffness matrix 
factorization [UTKU, MELOSH, SALAMA, CHANG (1986), UTKU, SALAMA, MELOSH 
(1986)]. The possible utilization of hypercube architecture is studied by NOUR-OMID, 
ORTIZ (1986) for the factorization and iterative solution of finite element equations. It 
is observed in these papers that the plain reduction of total number of operation does 
not always provide an optimal algorithm for vector and parallel processing. Increasing, 
number of sessions are devoted to discuss this aspect of finite element computations in 
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scientific meetings, and literature is being accumulated. However, further research and 
development still needs to be done before robust algorithms become available for a 
wide range of engineering applications which fully utilize the potential of modem 
computing machineries. 

Notable progress has been made in the last few years in the field of 
computational plasticity, in particular, design and analysis of integration algorithms for 
rate-independent plasticity constitutive equations. ORTIZ, SIMO (1986) summarizes 
the development of integration algorithms and investigates the accuracy in great detail. 
The basic idea is to generalize the return mapping algorithm based on the elastic 
predictor. An extension to viscoplastic constitutive equations are included in this 
paper. SIMO (1986) extends the results to finite deformation plasticity developing a 
simple and efficient finite element algorithm in which elastic deformation of finite 
amplitude is accommodated. 

Approaches to capture the localization due to the inelastic material response in 
global finite element computation have been developed which enable one to capture 
discontinuities occurring at subelement scale. The work by HUGHES, SHAKIB (1986) 
extends the conventional return mapping algorithm for /2-flow theory to the yield 

surfaces with comers. This extension can be used to capture the effects of localized 
plastic flow in a global manner. ORTIZ, LEROY (1986) proposes an algorithm which 
allows shear bands to generate inside of a finite element. The numerical results 
indicate the potential of the method to capture local failure of structures without 
excessive mesh refinement. Ideas to bring in the local failure mode in finite elements 
to capture localization effects under strain softening are discussed by WILLAM, SOBH, 
STURE (1986). A survey of the finite element technology for capturing the localized 
failure modes is given by NEEDLEMAN (1986). This is a relatively new field and 
further work will lead us to effectively calculate the localized microscopic material 
responses due to the embedded singularities in large scale structural systems. 

In the theoretical aspects of finite element development, the utilization of mixed 
variational principle has been and will be the main thrust of the research and 
development. This approach not only provides us with tools to improve the 
performance of elements but also lets us get a deep insight on the solution algorithms 
for linear and nonlinear problems. The mixed iterative methodology developed under 
HOST contract takes full advantage of these modern theoretical developments. Also, 
new ideas demonstrated in the MHOST program has attracted academic research 
interest in the computational mechanics community. 

In the application and computational aspects of finite elements, significant 
progress made in the last few years has not been fully incorporated in the MHOST 
computer code development. The information available in literature indicates that 
further development to utilize modern algorithms and coding technology would 
enhance the performance of the MHOST code considerably. 
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2. VARIATIONAL FORMULATION AND COMPUTATIONAL PROCEDURES 


Introductory Remarks 

This chapter is devoted to the problem statement and its solution strategy of a class 
of problems dealt with in the framework of MHOST finite element program 
package. We first establish the notation used throughout this report in the form of 
the differential equations and then derive the variational statement suitable for the 
construction of the mixed iterative finite element analysis. In order to maintain the 
generality, the nonlinear transient problems used as the vehicle for this develop- 
ment in which the conventional quasi-static problems are included as a subclass. 


The Global Solution Strategy 

The nonlinear problem to describe the inelastic response of a material body in an 
open, bounded domain ft with sufficiently smooth boundary dft is stated in this 
section and an augmented form of the generalized nonlinear Hu-Washizu 
variational principle is derived in the infinitesimal deformation setting. A generic 
solution algorithm is constructed which is valid for the quasi-static and the 
dynamic-transient analysis. 

The procedure discussed here can be incorporated with various solution 
algorithms and time integration operators other than the conventional Newton- 
Raphson and Newmark-13 methods. 

The equilibrium equation is written in terns of Cartesian components of 
stress tensor as: 

in a (21) 

where p is the material density, assumed constant, with a, a and / being the stress 
tensor, the acceleration vector and the body force vector, respectively. 

We assume a rate constitutive equation in a generic form 

u ( 2 - 2 ) 
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with D and £ being the material modulus and strain tensor respectively. The strain 
components are given by 


e »> *2 («iv + “>.*•)■*■ 4 


( 2 . 3 ) 


with being the initial strain due to the thermal expansion and creep effects. 
Usual equality constraints for the displacement components are imposed on the 
boundary such that 

Uj = iij on dQ! (2.4) 

and 

t J = o i jn l = tj on dQ 2 (2.5) 

The weak variational form associated with the above problem statement is 


i a u >"ij) = (p (fj - a j) •">) + (o° •“/) 


and 


( £*• ,\ojj - D ljkl {z k i - £*/)]) = 0 
e ij ~ 2 i u iJ + ® 


* 

a ijy 


( 2 . 6 ) 

( 2 . 7 ) 

( 2 . 8 ) 


where ( . , . ) denotes the usual L2 inner product over the domain, defined by 



uvdx 


ci 


m,v£^( 0 ) 


and <.,.> is the integral defined on the boundary defined by 

( u, v) = ^uvdx M ,v€L 2 (dQ) 

an 

with the superscript * indicating the arbitrary variations from an appropriate space 
of admissible functions. 
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Elimination of the stress and the strain from the above system of variational 
equations results in the virtual work equation in terms of displacement 

a(u,u*)~ (f,u*) = 0 ( 2 . 9 ) 



with a(.,.) being the usual energy product defined by 

a(u,v)=ju i jD i j k iV jJ dx 

a 

The essential boundary condition, Equation (2.5) can be incorporated by virtue of 
penalty 

a(o,n*) + e"'(« -!< 0 ,u*)= 0 (2-10) 

or equivalently 




( 2 . 11 ) 


It is obvious that the simultaneous weak variational statements. Equation (2.6 - 2.8), 
can be derived directly from the Hu- Washizu principle. This observation implies 
that the boundary conditions, Equations (2.4) and (2.5), are entered to the system of 
equations only via the conservation law for the linear momentum. In this setting 
imposition of boundary conditions for stress and strain is unnecessary. If such 
conditions are applied, the well-posedness of the problem may be disturbed yielding 
erroneous solution or possibly a rank deficient system of equations. 

Note that the penalization for the Dirichlet boundary condition does not 
require the space of admissible variation to fulfill die homogeneous counterpart of 
the same condition. 

Using the equal order interpolation function for all the variables involved 
in the analysis, we have 

«J = N k V jk (2.12A) 


a h j - N K a }K - N K u jK 


(2.12B) 
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4 - 

N k E,j K 

(2.12C) 

4 ~ 

NkS^k 

(2.12D) 

and for the 

data 


II 

N k Fjk 

(2.12E) 

g 0 * = 
fc */ 

- N k eI k 

(2.12F) 


resulting in a system of algebraic equations 


0 

0 

B 

f ^ 

u 


F-MA 



0 

D 

-c 

E 

> = < 

G 

► 

(2.13) 

B 

-c* 

0 

.s. 


10 j 




where the submatrices are defined as 


V=\V, L ] E=[E lkL \ S-Ivl 
E t j tKL = jN K jN L dx 

Cl 

C = [ O/JWKz] CijklKL = ^ N K N L dX 

Cl 

D = [ DjjkiKll A jklKL ~ K E l} kiN L dx 

Cl 

M = \\i ijia\ M ljKL = pj N K N L dx 

ci 

with the loading vectors defined by 

/ 7 = [/VJ FiL = $H L pf l dx + fNjjds 

n en 

G = [ Gjj J Giji = E, jkl e ° k! dx 

ci 

Note here, for the sake of simplicity, the integrated 


form of the rate constitutive 
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equation is assumed. We shall discuss the detail of incremental process and the 
stress recovery later in this chapter. The elimination of the nodal strain and stress 
lead to a displacement solution form: 


j(fi(C)' , Z)(C)' 1 fi') + P\V = F- Ma 


(2.14) 


with 

F= F+ {B(Cy'D(C)~ l G} (2- 15 ) 

whereas the standard finite element form based on the statement (2.9) is 
approximated by a somewhat simpler form 


{K + P)U = F-Ma 


(2.16) 


An iterative solution algorithm is constructed first to the quasi-static counterpart of 
Equation (2.13) assuming the time derivatives of displacement vector are 
sufficiently small: 

(i) Initialize the residual and displacement vector 

R:=0 ; 17:- 0 

(ii) Solve the preconditioning equations to update the displacement such that 


U:=U+ A~ 1 (f-R) 

(iii) Recover the nodal strain 

E:= C~ ] BU- £° 

(iv) Integrate the constitutive equation 

S: = jD{S,E)EdT 

T 

with t being the quasi-time associated with deformation history. 

(v) Evaluate the residual 

R:= F - BS ( 2 - 2 °) 

(vi) If the residual is small enough then exit; or else repeat from the step (ii). 


(2.17) 


(2.18) 


(2.19) 
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As is obvious from the above discussion, the choice of the preconditioner A is 
crucial to obtain the convergence characteristics necessary for the practical 
implementations. For instance, if we can construct 

A = {(fl(C)''Z>(Cf , fl') + P\ < 2 - 21 ) 

then no iteration is needed. However, the spareness of the finite element system of 
matrix is no longer exploited if the above form were to be utilized. As a reasonable 
compromise, we use the augmented equilibrium equation by the displacement 
stiffness matrix which is to set 

A = K + P C 2 - 22 ) 

Other methods of preconditioning has been investigated but so far no robust scheme 
is known for a wide range of solid and structural analysis. See for instance, 
MULLER, HUGHES (1984). Except for minor modifications, the present solution 
utilizes the form defined by Equation (2.22). 


Global-Local Analysis by Subelement Iterations 

Computational fracture mechanics aspects are discussed in this note of the inelastic 
analysis of turbine engine hot section components. The motivation of this 
endeavor is to seek an economically feasible numerical process without sacrificing 
the accuracy of the solution. 

The standard finite element method is to take the effect of embedded 
singularities into account by refining the mesh subdivision in the neighborhood of 
such points or alternatively to introduce special elements with singular functions in 
the same region. This currently available approach is often prohibitively expensive 
especially when it is applied to the analysis of a structure with multiple singularities 
of which each needs special treatment. 

The use of nonstandard schemes, in particular the mixed finite element form, 
is found in the previous section far more advantageous compared with the 
conventional schemes with application to regular nonlinear structural analysis. 
Hence its utilization for the problems including singularities is worth investigating 
to realize a major breakthrough in the computational fracture mechanics. 

The main objective to use the mixed finite element method for problems 
with embedded singularities is that higher resolution for stresses and strains is to be 
obtained without excessive mesh refinement for the displacement variables near 
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the singularities. 

In the framework of displacement finite element method, only available 
information primarily from the computation is the nodal displacements from 
which the strain /stress components are calculated by differentiating the shape 
functions. Then these quantities are used to evaluate the fracture mechanics related 
quantities such as the stress concentration factor and the J- integrals. 

Noting that the accuracy of strain /stress components is a full one order less 
than that of displacement for the displacement finite element method even without 
a singularity. Moreover, the approximation of strains /stresses could be extremely 
unstable near singularities and the area where the deformation is highly con- 
centrated. To obtain accurate results for a problem with embedded singularities by 
the displacement technology, the mesh refinement at around the singularities is 
unavoidable so as to obtain accurate and stable displacement field which indeed is 
the only source to generate the strain /stress approximations. For the best possible 
results, monitoring a posteriori error indicator and several passes of (adaptive) 
mesh refinement are advisable for such calculations. 

When the strains or stresses, sometimes both, are included explicitly in the 

finite element system of equations, in particular interpolated by cP- continuous 
basis functions, the resulting approximations are not only accurate but also stable 
having no spurious modes in the strain /stress recovery operator. 

For regular problems, the numerical representation of stress field in 
particular is quantitatively accurate due to the equilibrium condition being explicitly 
evaluated by the approximate stress field itself. Compared with the displacement 
method, a full one order improvement is indeed expected of the convergence rate 
for the mixed method for the same mesh subdivision. It is noticed that the role of 
displacement solution in the mixed form is to generate qualitatively the defor- 
mation mode from which indirectly the deformation gradient is extracted and fed 
into the stress recovery and equilibrium iteration operations. Therefore the quanti- 
tative measure for error in displacement vector contributes relatively less signifi- 
cantly to the overall evaluation of fracture related quantities. 

The approximate solution procedure for a deformable body with embedded 
singularities is first to solve the structural problem without crack which shall be 
specifically referred to as the global system and then to compute the deformation 
and the stresses near the singularity separately. This portion of the domain is called 
the local system. The concept is somewhat similar to the substructuring in the 
conventional finite element computations. 

Topologically the global mesh is constructed such that the minimum element 
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size is much larger than the scale of singularities so that no singularity is identified 
by the global system. The local system is constructed in a single globed element and 
coupled with the global system through common nodes between the global and 
local elements within. 

A major difference of the present approach from the substructuring is that the 
local refinement takes place a posteriori and the information extracted from the 
locally subdivided element is brought back to the global mesh via the residual load 
correction. This implies that no special algebraic treatment is required even when 
the material nonlinearity occur in the sub-element region. 

The solution strategy is indeed a mixed version of the hierarchic finite ele- 
ment process which has been investigated in the recent few years and the results ob- 
tained are encouraging for elastic stress analyses with singularities. The compu- 
tational strategy we propose here is to define the local mesh and the refined inter- 
polation altogether in the region near the singularity somewhat similar to com- 
bined h-p refinement. 

As the result, we shall obtain accurate stress-strain solution in the global 
approximation subspace enriched by the local approximations. 

The final solution to be obtained is in the global-local system with the global 
system being used only for preconditioning purposes. This means that the factor- 
ization of approximate global stiffness matrix obtained by the conventional displace- 
ment method may be sufficient to kick off the present solution procedure. 

Consider as the second example a single displacement element with an em- 
bedded singularity, a hole. Some computational aspects need to be considered here 
with application of present method to practical problems of turbine engine hot 
section components. From the programming point of view, a parametric repre- 
sentation of an elliptic hole by its size, orientation and location in a displacement 
element is convenient allowing the code to generate subelement mesh data when- 
ever this calculation is performed. The data structure could be simple as no perma- 
nent storage allocation is required for the subelement mesh. In Figure 2.1, an 
example of such a ready made subelement mesh is presented for a circular hole 
located at the center of displacement element mesh. 

On the other hand, it may be user friendly and perhaps conceptually more 
general to explicitly define the subelement mesh as an additional data set. The data 
structure then needs to be reviewed so as to allocate additional memory for sub- 
element mesh storage. 

As experienced in the h-version of adaptive mesh refinement, it is anti- 
cipated that the concept of subelement mesh could be used recursively in a similar 


S.Nakazawa 


March, 1991 


MHOST Theoretical Manual 


22 


fashion as the multi-level substructure technique. A well organized data storage 
scheme needs to be developed to realize the fully flexible implementation of the 
proposed scheme. In the following discussions, the subelement mesh is defined in 
the isoparametric element coordinate system rather than the physical space in 
which the global mesh is defined. 
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The stresses and strains are represented at the nodes of subelement defined in 
the first reference plane in an exactly same manner as the one dimensional 
example. 

We introduce variables to characterize the internal deformation of stress and 
strain subelement so as to obtain these quantity uniquely in an accurate manner. As 
mentioned in the first section, some further mathematical investigation is required 
to come up with an optimal combination of subelement mesh definition, functions 
to represent the deformation of subelement and the stress-strain interpolations. 

A major difficulty encountered in the multidimensional computation is to 
integrate the coupling matrix B over the subelement. Noting that once the size of 
hole and its location are identified, the values of displacement interpolation 
functions are obtainable at every nodal point of subelements. Denoting these 
quantities by N , we introduce the local displacement-type interpolation of N in 

these elements by 

N = N S k N k (2.23) 

and the coupling matrix is integrated approximately by 





(2.24) 


where the superscript S denotes the functions defined at the subelement level. The 

calculation is carried out on each local subelement and summed over the global 
displacement element. This general integration procedure plays a central role to 
couple the global and local mesh representing the residual load vector at the global 
displacement node entries yet reflecting the information at the local element level. 
For the consistent definition of this operator matrix, the transformation of shape 
function from global to local needs to be carried out as defined by Equation (2.23) on 
the isoparametric element space which is again mapped into another isoparametric 
space defined on each subelement. It is obvious that the location of nodal points 
needs to be given in the element coordinate system. 

When the location of subelement nodal points are specified in the physical 
plane, the inverse of nonlinear isoparametric mapping needs to be calculated to 
locate them in the element coordinate system. 

We outline the solution strategy in a general format. First we introduce the 
conventional stiffness equation with respect to the local hierarchial displacement 
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degree of freedom U as the preconditioner, which is 

K S U S =F S ; F S = F G -B S S G (2.25) 

derived directly from the virtual work principle in terms of displacement. A 
modified recursive form of mixed finite element equations is written again in terms 
of local unknown variables 


K S 

0 

B S 

U s 


F S + K S U S 

0 

D 

-c 5 

<E S 

> = < 

0 

B st 

-(?' 

0 

<*> 

> 


0 


(2.26) 


Setting tf (O) =0 and S °^ = 0 with the superscript being used for the iteration 
counter, we solve the system of displacement equation 


= t/OOS + (k 5 )' 1 ) f s -b s & n)S ) 


(2.27) 


Then the stresses and strains are updated in the subelement region by solving 
iteratively with u being the enriched displacement in the subelement region. To 
obtain this quantity, the factorization of displacement type stiffness equation is 
indeed unavoidable in this region. 

In the solution of local system of equations, the traction boundary condition, 
Equation (2.5) is specifically imposed which cannot be applied to the solution of 
global finite element process. In the part of global mesh where no singularity is 
embedded, conventional recovery procedure is used prior to the above calculation. 

Solution of Equation (2.27) is improved in comparison with the global mixed 
solution by the local information contained in the righthand side. The use of 
additional terms characterizing the deformation near the singularity needs to be 
considered. The simplest approach is to take the displacement as a dummy variable 
at stress-strain nodal points interpolated by conventional polynomial basis 
functions in a similar manner as the substructuring technique. 

When such variables are introduced the solution procedure becomes very 
close to what is called the multi-grid method in the finite difference context. The 
overall solution flow is: 

(i) Solve the global stiffness equations for the global displacement, i.e.. 
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lf:«tf + K G " 1 (F-.R) 

(ii) Recover the strain at nodes in the global finite element mesh. 

(iii) Evaluate the stress. 

(iv) Enter the inner loop and calculate the displacement quantity (used only as 
dummy variable) in the subelement mesh. 

U^DS.uOtS + ^-'jfS-BSsCW} 

(v) Evaluate the subelement stress and strain. 

(vi) Check the convergence in terms of residual in the subelement mesh. If not 
converged, repeat from step (iv). 

(vii) Evaluate the global residual including the stresses in the subelement. If 
convergent, start a new increment, or otherwise repeat from step (i). 

R:= F - BS S 

Through the inner loop, the global finite element mesh interact with the 
subelement and an accurate, high-resolution stress field is obtainable without 
increasing the size of global stiffness equations. 

The additional coding is required to incorporate steps (iv) to (vi), which is 
performed by adding a element routine to handle the embedded singularities, and to 
allocate the core storage associated with the arrays used for the inner iteration. The 
size of algebraic equations used for representing the embedded singularity is limited 
by restrictiong the mesh pattern and the level of mesh refinment introduced 
through the subelement technique. 
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Incremental Iterative Solution Algorithms 

The inelastic problem is solved through the deformation history in an incremental 
manner which is for a given state of displacement, strain and stress which satisfies 
the nonlinear algebraic equation (2.13) for a given load and displacement boundary 
conditions by UJE, S and F respectively. 


Adaptive Load Incrementation Procedure: 

For a given load increment A F, the increments of displacement, strain and stress 
denoted by A U, A E, and AS are determined in an iterative fashion as described in 
the previous section. Note in the preconditioning process the tangent stiffness K is 
used instead of the total stiffness equation. To follow a complicated equilibrium an 
automatic adjustment procedure procedure to control the size of load increment in 
die iterative process is introduced. The algorithm is schematically given as: 

(i) Set the residual vector R:= 0; 

Initialize the incremental displacement vector A U: = 0; 

Set the iteration counter /:= 0. 

Apply the proportional load A F:= XF with X being the current load factor. 

(ii) Update the displacement 

AU:=AU + dU (2.28) 

with the update 

dU:= K'\F 0 + AF-R) (2.29) 

being stored in the memory, where Fq is the load vector at the beginning of 
the current increment. 

(iii) For a given arc-length, find the total load factor update 

A: = A + dX (2.30) 

and the incremental load factor 

AA:=AA + dA (2.31) 

based on the spherical path formulation [Crisfield (1980)]. 
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(iv) Modify the load 

AF:= AAF (2.32) 

(v) Form the residual in the mixed manner as indicated in steps (iii) to (v) in the 
previous flow chart. If convergent, then exit; otherwise repeat from step (ii). 

BFGS Update Procedure 

For the iterative process using the mixed residual vector, any iterative method to 
improve the convergence characteristics applicable to nonlinear finite element 
calculations can be used for the present approximation method. For instance, the 
BFGS update procedure [Matthies, Strang (1979)] is utilized in the following fashion: 

(i) Initialize the residual vector R and the incremental displacement A 1 1 such 


that 

R:= 0 ; AU: = 0 (2.33) 

(ii) Modify the residual with i being the iteration counter 

(2-34) 

(iii) Intermediate displacement update 

dll = K~ l F (2.35) 

(iv) Complete the displacement update 

(2.36) 

(v) Form the new incremental displacement 

AU:=AU + dU (2.37) 


(vi) Form the new residual with respect to the updated incremental displacement 
using steps (iii) to (v) in the first flow chart. 

(vi) Check the convergence and if necessary repeat from the step (ii) or else, exit. 
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Note in the above algorithm, vectors v and w represents the iterative change in 
the residual and the displacement vector respectively so as to form the inverse 
BFGS update 

K~n = (l +w n + v„ wl) (2.38) 

It is possible to introduce the combined BFGS - arc length algorithms in the mixed 
iterative solution algorithms incorporating the line search technique. 


The Line Search Algorithm 

The line search algorithm is used for the acceleration of niterative solution 
processes. After an update vector du for the velocity is obtained for a current 
displacement increment, we seek for a search distance s which approximately 
satisfies the orthogonality condition 

g =5 u r R <fi - sb u) = 0 ( 2 . 39 ) 

An iterative algorithm [Abdel Rahman (1982)] is implemented to find zero for the 
above defined function g of which the extrapolation procedure performed at each 
iteration is illustrated in Fig.2.2. 



Figure 2.2 Linear Extrapolation in Nonlinear Line sraech 
The compuational procedure is schematically written as: 
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(i) Set S= 1 

(ii) Recover the residual vector 
R : = R (fi - sb u) 

(iii) Inner product to evaluate the function 
g =bu r R 

(iv) If g - e a predefined tolerance then exit; else: 

(v) Extrapolate linearly the search distance such that 


(vi) If s> Sg , with Sq being the prescribed upper bound, set s = . 


( 2 . 40 ) 


( 2 . 41 ) 


( 2 . 42 ) 
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Formulation of Finite Deformation Analysis Procedure with Particular Emphasis on 
Plasticity 

The formulation and algorithm for the analysis of finite deformation problem 
implemented in the MHOST mixed iterative method is described in this subsection. 
The formulation utilizes the kinematics of updated Lagrangian concept. The choice 
of deformation and stress measures in the finite deformation analysis depends 
primarily on the form of the constitutive description involved, along with the 
necessity that the measures are objective for rigid body motions. For elastic-plastic 
behavior where the response depends primarily on the current state of stress, a 
Cauchy stress and rate-of- deformation constitutive formulation is frequently the 
most convenient choice. 

The constitutive integration algorithms in the updated Lagrangian fromu- 
lation is essentially the same as those used for infinitesimal deformation cal- 
culations, with the exception that all tensors and integrals are evaluated with respect 
to the current configuration. The present version of this update strategy utilizes the 
nodal coordinates updated continuously during the iterations in conjunction with 
the deformation gradient evaluated simultaneoulsy at nodes. 

We shall now summarize the iterative solution algorithm for an increment. 
Following initialization, the process is repeated iteratively until the equilibrium has 
been satisfied. 


(i) 


Initialization: 


u' = IT" 1 S'=S M R = F 


(ii) Equation Formulation and Displacement Solution: 

A17' = K~ 1 r! 


(iii) Update Geometry: 

i / 1 * 1 = u' + \ir 


TrP-I + V 


(\ 




~A U' 
2 


(2.43) 

(2.44) 

(2.45a) 

(2.45b) 

(2.45c) 
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^J +1 = J + V( AIT> 1+ *- 

(2.45d) 


j 1+1 =^j + V 

(2.45e) 


^ I+ i = 1 

(2.450 


^ 1+1 =^J + V 

(2.45g) 


^ f+1 . J f+1 ( ti i+1 ) 1 

(2.45h) 

(iv) 

Strain Projection: 



A£ , ' +1 -(C , ' +1 )’ r B , ’ +1 AU i 

(2.46) 

(V) 

Stress Recovery: 



S' +1 = Or ^ ,+ ') T AE ,+1 ^' + ' + S' 

( 2.47) 

(Vi) 

Rotate Back to Global System: For stress components 



s M 1 (s. i **) r 

(2.48a) 


and for plastic strain components, we have 



(Epj i * 1 -* M MTV* , ) r 

(2.48b) 


Material tangent is updated as 




(2.48c) 


Internal variables susch as the back stress is also updated by 

(2.48d) 


S.Nakazawa 


March, 1991 


MHOST Theoretical Manual 


32 


(vii) Form Residual: 

K ,+1 = f(N i+1 f f i+1 N i+1 dx hl - B i+1 S i+1 (2.49) 

n i+1 

If converged exit, else update displacement via line search or displacement solution 
and repeat. In the finite displacement calculations, the stiffness array is assembled 
in the following manner and repeatedly updated in the iteration loop unless other 
iteration process (modified-, quasi-, or secant-Newton option) is specified. Note that 
this stiffness matrix update does not correspond to the classical Newton method in 

the mixed iterative solution framework because of the independent C^- continuous 
stress interpolation functions with K being the index for nodes, 

o h =N K S K (2.50) 

Whenever the stiffness matrix is reformulated in an updated Lagrangian 
algorithm, the shape functions, derivatives of the shape functions and determinants 
of Jacobians are evaluated with the nodal coordinates of the current configuration. 
In addition, the material and stress matrices are also the values computed in the 
immediately preceding iteration. 


Semidiscrete Finite Element Equations and Temporal Discretization 

The finite element equations of dynamic equilibrium for a deformable body is 
written in terms of nodal acceleration vector a, nodal velocity vector V and nodal 
stress vector S as 

Ma+ CV + BS = F (2.51) 

where M is the mass matrix defined by 

M = p^N T Ndx (2.52) 

Cl 

with r and N being the material density and the vector of interpolation functions, 
respectively. The damping matrix C may be defined as a linear combination of mass 
and stiffness matrices, M and K, such that 

C=piM+p 2 K (2.53) 
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with pi and p2 being preassigned constants representing the effects of viscous and 
structural damping, respectively. Symbolically, the stiffness matrix can be written by 




J(VN) r DVM* 

n 


(2.54) 


with E being the gradient operator and D the material modulus. The third term in 
Equation (2.47) represents the internal energy with B matrix being the discrete 
gradient operator defined in the finite element approximation subspace such that 



n 


Note that this array differs from usual displacement-strain matrix in conventional 
finite element computations. 

Here the equal order interpolation using the same basis function is assumed 
for the acceleration, velocity, displacement, and stress. Introducing the process to 
recover the nodally interpolated strain by 


E = GBU 


(2.56) 


where G is the diagonalized gramm matrix of which Ith diagonal entry is calculated 


by 


G = 


Ijw 


Ndx 


k- in 


(2.57) 


with N no d e being the total number of nodes in a mesh. Then the stress is 
recovered at nodes. 


Note that all the integration appeared in the above equations are evaluated 
approximately using the numerical quadrature. The particular choice of the 
numerical integration rule is discussed in the subsection of element formulation. 

Now we construct an abstract iterative process for the semi-discrete Equations 
(2.52) using the displacement preconditioning. The approximation of internal 
energy by the product of stiffness matrix and total displacement vector u leads to a 
recursive expression 

mJ 11 + CV '* 1 + KU M = F- (BS'-KU') (258) 
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with subscripts denoting the iteration counter. For the linear elastic problem, the 
displacement converges linearly for an appropriate initial condition such a s U = 0. 

Note that the strain projection procedure and stress recovery operation, 
Equations (2.56) and (2.58), are executed every time when the displacement vector is 
updated. 

The fully discretized system of equations is derived for the mixed iterative 
solution by introducing the temporal approximation. For the sake of simplicity and 
clarity when it is implemented, the Newmark family of algorithms is introduced in 
a finite difference fashion e.g., HUGHES (1984) such that the dynamic equilibrium is 
satisfied at time t + At 


M Of*# + CV t +£i + BSt+&t ~ F t +£ (2.59) 

which leads to a time discrete displacement preconditioning 

= F -(BS t +u l - KU t+ x) (2.60) 


and the updates for displacement and velocity are given by 


2 

I/,'.* - U, + MV, + ^-{(1- 2 PV, + 21^4 


(2.61a) 


v; +A , = V t + At j(l - Y )a t + ya \ +dl ) 

(2.61b) 

Introducing 


A U\ = UUm ~ U t 

(2.62a) 

and 



5i/; = au; +1 - au\ 

(2.62b) 

we obtain a recursive form 


< 

1 V 'I 1 j 

— i-r M + -77 C + K\bU 
PAf 2 J 
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BSl 


t+Ai 
f 




{ ( ^ A 


llPAf- 


Af + 


' Y ' 


PAf 


c Ai/; 


^-ilcv+ 


X-1 

2(3 . 


AfCdj + 




( 


MV t + 




Mat (2-63) 


with the incremental displacement Du being updated by 

ai/; = au;+5i/; 

and the incremented strain recovered at node by 


(2.64) 


AEi.a = G _1 B T AUU n (2-«) 

In a general setting of history dependent inelastic constitutive integration, the stress 
recovery process is written as 



( 2 . 66 ) 


which is fed into the equilibrium iteration defined by Equation (2.63). 


Computational effort involved in the above iterative process for direct time 
integration of the discrete equation of motion is the assembly of element operator 
matrices appear in the left-hand side of Equation (2.63). This is exactly the same 
amount of effort required to solve the problem by the conventional displacement 
method. Additional computation required for the mixed solution is a few back 
substitution and strain and stress recovery at nodes which are relatively 
inexpensive compared with the matrix assembly and factorization. This additional 
computational cost would vanish when the scheme is applied to nonlinear 
problems. The matrix assembly and factorization or back substitution needs to be 
performed repeatedly no matter which finite element method is used to drive the 
solution. Therefore, seemingly complicated derivation of the mixed iterative solu- 
tion scheme does not increase the computational effort in comparison with the 
displacement formulation to fully enjoy the advantages of equal order interpolation 
of displacement, strain and stress. 

With particular application to turbo machinery blades and other rotating 
structures, the stiffness matrix is modified to include die effects of stress stiffening 
and centrifugal mass terms. These are linearized models of large displacement and 
follower force around a given state of stress due to the linear elastic response of 
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structures under a given angular velocity. 

The centrifugal load vector is generated by integrating the body force 

io = J*Np co 2 r(x) dx (2.67) 

n 

with (0 being the angular velocity and r the distance between the point and the axis 
of rotation. Under the centrifugal loading given above, the mixed iterative 
algorithm results in a certain stress field represented by the nodal stress vector s 
such that the discrete equilibrium equation is satisfied. Our interest is to calculate 

the linearized response of structure under the state of stress given by S®. Taking 
the initial stress terms into account, the equilibrium equation is modified (excluding 
the damping term) to 

Ma+ K^U + B t S = F (2.68) 

where K is the geometric stiffness term associated with the initial stress field s 
given by 

Kq= ^N T SVNdx (2.69) 

n 

Note that Equation (2.64) is a general expression for the linearized response of 
structural systems under a prescribed initial stress field. 

When this term is included in the mixed iterative computations, the abstract 
recursive form. Equation (2.55) is modified to 

Ma+ (Kc + K)U m - F - {B T S i - (xg + K)U,) (2.70) 

or in terms of displacement update 

(Kg + K)l/ M =F-(B T Sj-KgU , 1 ) (2-71) 


with 


If ,♦! = U, + AUj 


In the residual force calculations, the initial stress matrix needs to be evaluated 
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repeatedly with the stress field calculated at the beginning of the increment. 

The additional force due to the relative motion from the equilibrium 
position referred to as the centrifugal mass term is also added to the equilibrium 
equation resulting in 

Ma+(K%-M C )U M + B T S = F (2.72) 


where M is the centrifugal mass matrix defined by 



u 2 N T LNdx 


(2.73) 


with L being the matrix of unit vector m parallel to the axis of rotation given by 


1 = 


2 2 

m 2 + m 3 

-m 1 m 2 

~m 3 m l 


-m 3 m 2 

2 2 

m 3 + mi 

-m 2 m 3 


-m 3 m 1 

-m 2 m 3 

2 2 

mi + m 2 


(2.74) 


Note that when this term is activated in the quasi-static and transient dynamic com- 
putations, the same correction of the residual vector calculation as for the stress 
stiffening terms as indicated in the right-hand side of Equation (2.71). 

Typically, these optional terms are activated for the vibration analysis of 
rotating structures by extracting the modes from an eigenvalue problem 

j M - x(Kg - M c + K)| X = 0 (2.75) 

In the modal analysis, the displacement preconditioner is assumed to 
represent the stiffness of structures. However, the stress data used in the evaluation 
of initial stress terms are calculated by the mixed iterative method which generates 
more accurate stress field than the conventional displacement method. Hence, the 
solution of Equation (2.75) is expected to be somewhat more accurate than the 
displacement method solution. 

For die detail of the derivation of initial stress and centrifugal mass matrices 
and fully worked out examples, see THOMAS, MOTA, SOARES (1973), RAWTANI, 
DOKANISHI (1974). 
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Additional Remarks 

In the present implementation of mixed iterative solution, the displacement 
stiffness equations need to be assembled and factorized. To enhance the efficiency of 
MHOST code, the profile solver replaces the constant bandwidth Crout decomposite 
algorithm. The implementation of solver is based on the code published by 
TAYLOR (1985). The code presented in the paper include both symmetric and non- 
symmetric cases using the Crout decomposition. The implementation of MHOST 
code excludes the nonsymmetric part of code to gain the maximum efficiency. The 
profile storage of the global stiffness equations reduces the memory requirement 
and time required for factorization. Computational overhead is added to prepair the 
elimination table for a stiffness matrix stored in profile form. 

In the quasi-static computations, option added to accelerate the rate of 
convergence for nonlinear iterations are made functional incorporating with the 
profile solver. These options are the modified Newton, quasi-Newton method of 
inverse BFGS update, secant Newton implementation of Davidon rank-one quasi- 
Newton update, conjugate gradient and line search. The spherical path version of 
arc length method incorporating with the modified Newton iteration is also made 
available for the profile solver. 

The eigenvalue extraction procedure implemented in the MHOST code is 
based on the subspace iteration documented by BATHE, WILSON (1976). In this 
process, a large eigenproblem of structural systems is mapped into a subspace of 
finite dimensions and the Jacobi iteration is performed to solve the small eigen- 
problem set in the subspaces. This eigenvalue extraction procedure is used to 
compute the natural frequency and vibration mode of unstressed and prestressed 
structures and the buckling analysis of prestressed structures undergoing elastic or 
inelastic deformation of infinitesimal or finite amplitude. 

Note here that in the eigenvalue extraction, the displacement stiffness matrix 
is used to represent the response of structural systems. However, in the buckling 
and model analyses with prestress, the improvement of stress field generally 
improves the quality of solution considerably. The utilization of mixed stress inter- 
polation in the eigenvalue analysis can be viewed as a perturbed eigenvalue 
problem in which the effect of independent stress approximation is regarded as the 
perturbation to the displacement finite element system. It is anticipated that an 
algorithm proposed by SIMO (1985) and implemented for the NESSUS probabilistic 
finite element code may be usable for eigenvalue extraction of mixed system of 
finite element equations. 
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3. ELEMENT TECHNOLOGY 


Introductory Remarks 

The MHOST program consists of a library of highly accurate linear quadrilateral and 
hexahedral elements. The selective reduced integration and assumed stress 
interpolation are used to improve the element responses under bending load. To 
take the full advantage of these formulations, a rational local coordinate 
transformation derived from the polar decomposition of the Jacobian matrix is 
developed and implemented in the code. This highly sophisticated algorithm 
makes the response of element insensitive to the isoparametric distortion. 

For the sampling and nodal projection of element strain components, the 
trapezoidal integration rule is utilized for the maximum accuracy and stability. The 
oscillation of strain /stress fields often observed at the element integration points are 
filtered out by this procedure. The resulting mixed iterative finite element method 
is capable of producing very accurate displacement, strain and stress simultaneously . 

The nodally continuous stress assumption used in the beams and shells 
improves the accuracy of MHOST solution cosiderably in comparison with the 
conventinal displacement strategies. For example, the nodally exact displacement 
and moment solution is obtained for a cantilever beam problem modelled by the 
MHOST shell element subjected to a point load in the transverse direction at the tip. 
This shows the capability of the mixed iterative shell formulation to capture the 
linear moment field exactly by using only the linear finite element basis functions. 

To avoid numerical instabilities commonly observed for the four node shell 
elements such as the numerical locking and hourglass displacement modes, the 
stiffness matrix for the shell element uses the selective integration for the 
transverse shear terms and a simple and efficient hourglass control scheme 
proposed by BELYTSCHKO, TSAY & LIU (1981). 

The sophistcated three dimensioned element implemented in the MHOST 
program exhibits superior accuracy in comparison with the conventional 
isoparamatric 8 node brick element derived directly from the displacement method 
formulation. The selective shear integration /assumed stress formulation in 
conjunction with the local coordinate trasformation based on the isoparametric 
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mapping results in a highly accurate displacement field even at the initial 
displacement preconditioning phase without iterative improvement. The iterative 
process drives the solution further toward the direction of high accuracy for the 
given mesh. 

To generate a good displacement update for preconditioning purposes, two 
major refinements are incorporated in the present computational procedure. These 
are the improved version of filtering scheme for the selectively reduced integration 
and the modified numerical quadrature for plates and shells to avoid possible 
kinematic mode excitation. 


Coordinate Transformation and Filtering Algorithms 

An element coordinate system is developed in this section which neutralize the 
rotation associated with the isoparametric representation of the element. The pro- 
cedure is similar to the kinematics of finite deformation processes. It is developed 
to handle two- and three-dimensional elements in a unified manner. 

The isoparametric transformation is defined by the Jacobian 



where x denotes the physical coordinate vector with X being the reference 
(isoparametric) coordinate system. The isoparametric transformation is 
decomposed to the rotation and stretch by virtue of the polar decomposition 

/ = RU 


where R denotes the rotation tensor with U being the right stretch tensor with 
respect to the isoparametric mapping. The rotation tensor is calculated by 

R = ]U~ l 

Using this orthogonal tensor, the rotation neutralized coordinate system X is 
defined by 


X=gx 


or, by using the indicial notation, 
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y m - o m y 

A ~gk x k 

where 

g-R ' 1 

The figure depicts the definition of coordinate systems used in this development. In 
order to maintain the flexibility in the development of numerical integration 
procedures for element arrays, the operation is carried out only at the element 
centroid. The coordinate system defined at the element centroid is used for the 
definition of vector and tensor components in the element. 

The right Cauchy deformation tensor is defined by a product of Jacobian as 

c = /7 

or alternatively, in terms of the right stretch tensor 

C=U T U 



a) The Parent Element in the b) A Typical Element in the Physical 

Isoparametric Space Plane 

Figure 3.1 Definition of Coordinate Systems for a Linear Isoparametric Element in 
Two Dimensions 


In terms of the right Cauchy deformation tensor, the right stretch tensor is expressed 
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as 

U~ 1 = C~' = Nc^N c 

Hence, the rotation tensor is calculated from the expression 

R = /( N t c \^N c ) 

Eigenvalues and eigenvectors of the right Cauchy-Green tensor are calculated by 
using the Cayley-Hamilton formula in two dimensions whereas an iterative process 
is used for the three dimensional application of this coordinate transformation. The 
Cayley-Hamilton formula used in the polar decomposition, is given as 

U = -j=L==(c + 

Jl c + 2V3c 

where 



I c = trac£,n c = detC, 


and 



0 

1 


The Jfc/-th component of the linear infinitesimal strain tensor in the local 

coordinate system e (W ] is written in terms of components in the global coordinate 
system such that 


(*!>„ JOJD.OJ 

ij 5 i Sj c 


Summation is not implied the superscripts in the parentheses. Using the strain 
components in the global coordinate system, we have the strain component defined 
in the local coordinate system by 


(k> = (Q (1) 
c &m 5rt °rm 


Hence the ij - th strain component can be represented by the product of coordinate 
transformation and the original global strain components by 
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un mu) «) <n 

5 , 5 j 5 m on mn 

It is interesting to note that, the full coordinate transfomation system is an identity 
tensor 

(«(/>(*>(/> = j 
Si Sj 5m Sn A 

Introducing the filtering matrix for the selective integration defined by 
r {k) _ „(*)„ (k) 

G , *j 

The strain component in the local coordinate system can be written in a simpler 
manner 

e un = G k G ‘ s 

ij im jn mn 

The volumetric strain in the element coordinate system is defined as 

6 <V) =S** 

- r k r 

inf* jn ^ 

Similarly the deviatoric strain components in the element coordinate system is 
written as 


e U) = gtt 

-G.VO-rf'X 


iw jn 


mn 


In the numerical integration process, at the full integration points we evaluate 

e = ae (v) + 138 (</) 

whereas at the reduced integration point, 

e = (1 -cOe <v) + (l- (>)e < ‘ ,) 
typically with 

a 3 1 ; P “ 0 
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This effectively filters the sprious bending stiffness associated with the linear iso- 
parametric element formulation based on the conventional displacement formu- 
lation. 


Defining the original strain component E by 


E = BU (3.15) 

with B being the usual strain-displacement matrix calculated at the quadrature 
points and U the nodal displacement vectors, the energy functional 1(U) is 
obtained in terms of nodal displacement for linear elasticity as an example by, 

I(l/)= j [f A E (v) ® E (v) dx + Jm E ( d) ®£ (4> <£r 
ta o 

(3.16) 

where <8> denotes tensor inner products, A and p are the Lamme-Navier coefficients 
and F and T are the prescribed body force and surface traction respectively. 

In the above energy functional, either the first term or the second term is 
under-integrated to realize the necessary effects of the selectively reduced integ- 
ration. 

In computations, it is often more convenient to construct the strain-displace- 
ment matrix selectively, so that the above energy functional can be written as 

1(U ) = | DB S dxU ~ fT U (3. 17) 

n 

where F is the collection of prescribed load terms which appeared in equation (2.48). 
The new strain- displacement matrix B is the selectively sampled matrix equivalent 
to a certain mixed method under the isopara me tric distortion. In the matrix 
notation, the filter for the element volumetric strain is written as 

£ {V) = G T IGBU (3.18) 

and for the element deviatoric strain 


-^F'Udx- ^T*Uds 

n ad 
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£ {D) = G t (1~ I)GBU 


(3.19) 


at each integration point, with 1 and I being matrices 



'l 

l 

l" 


'l 

0 

o' 

II 

l 

l 

l 

II 

o* 

II 

0 

1 

0 


.1 

l 

l. 


.0 

0 

1 . 


Hence, we have 

B s = G T IGB + G T (1~I)GB 


(3.20) 


(3.21) 


where either the first or the second term is sampled at the reduced integration 
points and substituted to the array associated with the full integration points. Such 
an operation is trivial only for the linear Lagrangian elements with a single reduced 
integration point. 


Hourglass Control Algorithm for Plates and Shells 

The kinematic model generated by the reduced integration of transverse shear terms 
in the 4 noded bilinear plates and shell elements often induce numerical noise in 
the displacement preconditioning operation. However, excitation of these modes 
often referred to as the hourglass mode does not significantly deteriorate the 
convergence of the mixed iterative algorithms. To eliminate these modes filtering 
methods have been developed either to constrain the modes by modifying the 
stiffness equation (a priori hourglass control) [BELYTSCHKO, FLANAGAN (1982)] 
or to filter out the sprious noise after the nodal displacement is obtained ( a 
posteriori hourglass control) [ODEN (1983)]. 

Noting that the fully integrated stiffness matrix does not contain the 
kinematic modes, we utilize a simple version of a priori hourglass control algo- 
rithm originally developed by BELYTSCHKO et.al. (1982). Note also that the fully 
integrated equation locks the structure which produces excessive strain energy 
associated with the transverse shear tenns for a given displacement field. 

Constructing the transverse shear stiffness matrix K by adding 

K s - e/4 2 * 21 + (1 ■ - e) K% xl) P- 22 ) 

where e is the parameter associated with the aspect ratio of the element, the 
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kinematic mode can be eliminated without the expense of numerical locking, due to 
the insignificant participation of the fully integrated terms. 

The formula to calculate the filtering participation of K is calculated at the 
centroid of the element such that 

e = e 0 t c /h (3.23) 

with t being the thickness of the element at the centroid and h is the mesh size 

fc=|(l*l-*4M*2-*3D < 3 - 24 ) 

Typically e 0 = 0.01 is used for the validation and verification exercises. 

No additional computational cost is involved in this hourglass control 
algorithm because the element stiffness equations are integrated selectively and K is 
readily available without adding a new integration procedure. 


Assumed Stress Element Formulations 

In order to boost the performance of continuum elements, a family of assumed 
stress element formulations is added in the mixed iterative solution procedure. 
Regular continuum elements often lack the appropriate deformation modes to 
model shell-like structures in a satisfactory way. This was observed by AHMAD, 
ZIENKIEWICZ (1972) in the development of the 8-node thick shell element. The 
problem was circumvented by a selective reduced integration formulation. There 
are a class of shell elements incorporating special interpolations for the transverse 
shear terms to retain accuracy for the bending behavior of shells, e.g.. Heterosis 
element of HUGHES, TEZDUYAR (1978) and the 8 and 9-node thick shell elements 
proposed by HINTON, HUANG(1986). Four and eight node versions of such 
elements are proposed by NAGTEGAAL, NAKZAWA, TATEISHI (1987). 

The original library of continuum elements implemented in MHOST 
(element types 3, 7, 10 and 11) are based on the selective reduced integration with the 
coordinate transformation in order to retain invariance with respect to the 
isoparametric tmsformation. In recent paper by PIAN, SUMIHARA (1986), a 
hybrid stress element which exhibits excellent bending behavior for distorted 
configurations. The element uses five independent stress parameters to define the 

element stress field. 

The assumed stress element formulation developed in this program is a 
family of continuum-type elements with enhanced bending behavior with 
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reasonable accuracy even at a high aspect ratio. The elements were constructed 
using tin assumed strain formulation. The basic strategy is:(i) the identification of a 
set of independent stress modes representing the desired element behavior, (ii) the 
construction of a corresponding set of strain modes which results in a proper set of 
stress modes. The formulation is the conventional strain-driven form and imple- 
mented in the existing MHOST framework for material nonlinear problems. The 
strain modes are used to interpolate the element strain field and are projected to the 
displacement gradients by a mean square manner. All assumed strain modes are 
expressed in terms of a local element cartesian coordinate system obtained by polar 
decomposition of the isoparametric mapping at the centroid of the element. The 
stretch tensor obtained in the polar decomposition is used for computing scale 
factors to make the element computation dimensionless. 

The library of assumed strain elements implemented in MHOST includes 4- 
node quadrilaterals for plane strain, plane stress and axi symmetric problems, and an 
8-node solid element for modeling three- dimensional continua. The current 
implementation can be used with either the displacement or the mixed formulation 
and supports different integration rules for strain projection and residual recovery. 

Additional Remarks 

LAn algorithm is adapted in MHOST to allow the nodal definition of pressure 
loading on 2 and 3D continuum meshes. The algorithm is based on a nodal 
assembly of tributary areas at each node in such a way that a unique outward 
boundary normal vector is defined at each surface node. These normals define the 
effective surface orientation and the direction of the applied pressure at the node. 

For small deformation problems, this operation is carried out only once, 
during the first element assembly loop, and the resulting boundary normals are 
used to compute consistent pressure loads throughout the analysis. 
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4. ITERATIVE SOLUTION ALGORITHMS 


Introductory Remarks 

In this section, iterative solution algorithms used for the solution of mixed finite 
element equations are discussed in deatail. The fundamental solution strategy is 
presented in the following flow chart. This process is simlar to the modifed Newton 
approach for the nomlinear finite element calculations. One can find simlarity of 
the process to the Uzawa's algorithm for constrained minimization problems. 

Step 1: Initialization 

r n = f and u n = 0 

Step 2: Nodal displacement update 

||rt + 1 = u” + K -1 T n 

Step 3: Strain projection to the nodes 

e /!+l =C -1 E u n+1 

Step 4: Stress recovery at the nodes 


cjfl+l = CT t G e* +1 


Step 5: Form the residual 


rfl+1 = f -E t o n+i 


If || r n+ l |, > tolerance go to Step 2. 

where notation defined in the first section is used to identify vectors and matrices 
appeared in the mathematical expressions. 

The displacement stiffness matrix K is used algorithmically to (a) obtain the 
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initial trial solution used to start the iterative procedure for the mixed problem, and 
(b) help iterate the intermediate results towards the "exact" solution for the mixed 
system of equations. 


The Davidon Rank-one Secant Newton Update Procedure 

The role of the displacement stiffness K in the solution of the mixed problem is 
essentially algorithmic. Therefore, we can use a simple update procedure to try to 
improve our matrix between iterations. One of the simpler procedures of this type 
is known as the Davidon rank-one secant Newton update. This procedure is equi- 
valent to replacing the K matrix in Step 2 of the basic algorithm by an iteratively 
updated matrix of the form 

^-1 _ v -\ A (du^-K^ 1 Y") ® (du^-K^ 1 Y") 

" 1 ( du” - K^ 1 y” ) • Y n 


where 


yn = r n - r n ~l 

i.e., the difference between the residual vectors obtained for two consecutive 
iterations. Of course, this is not how we do it in the program! As we shall see, 
there is a far more efficient way of programming this method. 

The best way of programming this algorithm will involve only a small 
number of additional vector operations when the displacement updates are 
computed. Using this approach, we will not even need to explicitly store the update 
vectors employed by the algorithm! An efficient implementation of the Davidon 
algorithm can be summarized as follows: 

Step 1: Initialization 

r n = f and du” _1 = du” _1 * — 0 

Step 2a: Solve for the trial displacement 

du n = K _1 r n 

Step 2b: Compute the following terms 
y n — fn — |*w— 1 
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_ ( du^ - ! * + du" - du" -1 ) • r" - l 
( du"- 1 * + du" - du"- 1 ) • y” 

b = -c 

a= 1 -c 

Step 2c: Update the nodal displacements 
du" * = a du" + b du" - l * + c du"~l 
u" + 1 = u" + du" * 

Step 3: Strain projection to the nodes 
efl+1 =C _1 E u" +1 
Step 4: Stress recovery at the nodes 

CT" + l=C _T Ge" +1 

Step 5: Form the residual 
r" + l =f-E T a" +1 
If ;J r" +1 || > tolerance go to Step 2. 


As we can see, the only difference from the basic algorithm is in Step 2a-c. Please 
notice that the first pass with the Davidon algorithm is essentially the same as in 
our basic algorithm. In a way, the update procedure is "turned on" by the second 
pass through the iteration procedure. This algorithm can also be used in 
conjunction with line searches. The procedure for combining the two algorithms is 
very straightforward, and we will not need to complicate our discussion with those 
details. 


Solution Algorithms for Eigenvalue Problems 

Eigenvalue problems arise in a number of different situations. Let's start with the 
very familiar dynamic eigenvalue problem for (small amplitude) undamped 
vibration of a linear elastic structure. This will give rise to a set of matrix equations 
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of the form 


(K-XM)({) = 0 

wher M and K are the consistent mass and displacement stiffness matrices for the 
problem. When the stiffening effects due to the initial stress field and/ or the 
centrifugal mass effects observed in rotating machinery are included, the eigenvalue 
problem will become 

( ( K + Kg - M w )-XM)<j)=0 


where 



B T 0°B</n 


^ = 0)2 pN T H T HN da 

Jn 

are the geometric stiffness and the centrifuged mass matrix. 

The first one of these matrices is also used to set up the buckling eigenvalue 
problem for a linear elastic structure. If the problem is linearized about the origin 
(the usual case), the matrix equations will be 

(K-XKg)<t>=0 

A related eigenvalue problem comes up in nonlinear analysis when we are trying to 
detect the presence of a nearby bifurcation. This problem can be expressed by a 
matrix equation of the form 

( K 77 - X Kg ) <j) = 0 

where 

K n = I B T D n BdQ + B t o n B dQ. + other possible terms 

Jn Jn 
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B t a n B da 


are some tangent material and geometric stiffness matrices, both linearized at the 
end of the last load increment. The presence of a bifurcation would be indicated by 
the vanishing of one or more eigenvalues. 

Finally, we must complete our list with the deformation modes problem which is 
an eigenvalue problem of the form 

( K-A I )4> =0 

where I is an identity matrix with the appropriate rank. This problem may be used 
to assess the quality of the finite element approximation for certain classes of 
problems. It hardly ever is used in actual engineering practice. 

The subspace iteration procedure is employed in the MHOST finite element 
program to solve the generalized eigenvalue problem. We will present this 
algorithm in terms of the dynamic eigenvalue problem although, in reality, the 
subspace iteration algorithm can be used to solve any of die eigenvalue problems 
discussed in the preceding pages. The basic implementation of the subspace 
iteration procedure will involve the following steps: 

Step 1: Select the initial trial vectors 

<!>«= t4>?4>2 <t>ml 

Step 2: Solve for the subspace vectors 

R” = K -1 M <t> n 


Step 3: Construct the reduced eigenproblem 

K” = R” t KR" 


M^R^MR" 

Step 4: Solve the reduced eigenproblem 

(K"-A” M n )ip? = 0 
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for i - 1, 2 ... m using the Jacobi iteration method. 
Step 5: Obtain the improved eigenvectors 

<t>" = R" [ \f% < ] 

If not converged go to Step 2 
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APPENDIX COMPUTER CODE ARCHITECUTRE 


Program Organization 

The architecture of the MHOST code is schematically shown in Figure A.l. The 
execution supervisor routine (SUBROUTINE HOST) controls multiple analysis 
modules in a consistent manner. In Version 4 of the MHOST code, a pair of 
subroutines are added to check the consistency of the Parameter Data and to generate 
internal flags for the selection of analysis modules. The intention is to avoid 
possible occurrence of users to combine die features of the program package not the 
way the developers designed. 

The analysis modules represents the control structure of incremental 
iterative algorithms implemented in the MHOST code without going into detail. 
The source program with comments is designed to serve as die schematic flow chart 
of the computational process. 

The schematic flow of the element and nodal data manipulation is coded in 
the element assembly submodules which are entered from analysis modules. The 
element assembly submodules performs operations independent of element types 
and constitutive models. The assumptions introduced in the mechanics aspects of 
the formulations are explicitly coded at this level. These submodules call the 
librarian routines for elements and constitutive models. 

The element librarian subprogram (SUBROUTINE DERIV) generates quantities 
unique to the element used in an analysis. The MHOST code uses the standard 
finite element matrix notation as in ZIENKIEWICZ(1977), ZIENKIEWICZ & 
TAYLOR (1990). The element specific information returned from the librarian sub- 
program is the strain-displacement array referred to as B matrix in the previous 
subsections. 

The constitutive equation librarian subprogram (SUBROUTINE STRESS and 
BMSTRS) is designed to incorporate the nodal storage of stress and strain values. 
The STRESS subroutine sets up the loop over the points at which constitutive 
equations are evaluated. The current implementation is a nested double loop with 
outer loop being over the node and inner loop over the integration layer through 
the thickness. Note that the conventional displacement method can be recovered by 
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restructuring this loop in conjunction with a few minor modification in the core 
allocation for stresses and strains. 


Execution Supervisor 


Level 1 


Analysis Driver 


Level 2 





Element Library 

Level 3 




Constitutive Equation Library 

Level 3 




Work 

Group 

Subpro 

grams 


Algebraic Operation Utilities 


Level 4 


Figure 1 Organization of MHOST Program 


The librarian subprogram calls the constitutive equation package (SUB- 
ROUTINE NODSTR) from which individual subprograms for initial strains, stress 
recovery and material tangent are called. The librarian subprograms also controls 
the pre-integration of stresses, strains and material tangent over thickness for the 
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shell element. 

Note that the evaluation of the constitutive equation is one of the most costly 
operation in the nonlinear finite element computations. An attempt is made to 
minimize the execution of this process during the incremental iterative analysis, in 
particular the recovery of the residual vector. At the beginning of each increment, 
this process may be executed to evaluate the material tangent and contribution of 
initial strain terms which are often necessary for proper displacement pre- 
conditioning. 

Except for a small amount of information related to the convergence of 
iterative solution, all die report generation is performed at the end of increment. 
Optionally, post-processing and restart files are written at the end of user-specified 
increments. Generic reporting subprograms are called from analysis modules inside 
of the loop over the increments. 

A brief description of major subprograms are given below: 

Execution Supervisor: 

MAIN PROGRAM - declares the work space in blank common as an integer 
array. Also defines system parameters which are machine independent. 
Then pass the control to actual execution supervisor SUBROUTINE HOST. 

SUBROUTINE HOST - controls the sequence of execution of the analysis 
modules. First, this routine executes the control parameter data reader 
SUBROUTINE DATIN1 and check the consistency by entering SUBROUTINE 
LETCMD and RUNCMD. The current structure of code allows certain 
combination of two analysis modules to be executed sequentially. Analysis 
modules called from the execution supervisor are: 

SUBROUTINE STATIC for quasi-static incremental iterative solution. 

SUBROUTINE DYNAMT for the transient tome integration of dynamic equili- 
brium equation in an incremental-iterative manner. 

SUBROUTINE MODAL for the eigenvalue extraction for vibration mode 
analysis. This subsystem may be executed after the quasi-static analysis for the 
modal analysis of prestressed structures. 

SUBROUTINE BUCKLE for the eigenvalue extraction for buckling load 
calculation. This subsystem is executable only after the quasi-static analysis. 

SUBROUTINE FRONTS for the quasi-static incremental iterative solution by 
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the out-of-core frontal solution. Note that certain options are not available in 
this subsystem. 

SUBROUTINE SUPER for the linear dynamic response calculation by the 
method of mode superposition. 

Input Data Reader : 

There are three major subprograms: 

SUBROUTINE DATIN1 - Reads and interprets the parameter data input called 
by the execution supervisor. All the default values for control variables are 
set in this routine. 

SUBROUTINE BULKIN - Reads and interprets the finite element model 
definition data and prints the mesh and loading data for the initial increment 
(number 0). The bulk data reader is entered before the execution of analysis 
modules. This subroutine enters following lower level routines: 

SUBROUTINE INITH for the memory allocation of integer workspace in the 
blank common to store nodal and element data. 

SUBROUTINE DATIN2 for the model data input. 

SUBROUTINE DATOU1 for reporting the model definition data. 

SUBROUTINE CHKELM for the detection of clockwise element connec- 
tivity definition which results in the negative Jacobian matrix. This 
routine automatically corrects the connectivity table to counterclockwise 
direction. 

SUBROUTINE SUBDIV for the memory allocation for the subelement mesh 
data and automatic mesh generation for the subelements. 

SUBROUTINE INCRIN - Reads and interprets the loading and constraint data 
for each increment. This routine is invoked by individual analysis modules 
from the inside of loop over the increments. 

SUBROUTINE DATIN3, a small subset of the bulk data reader SUBROUTINE 
DATIN1 is used for actual operations including the initialization of arrays 
at the beginning of increment. 

Algebraic Operation Subsystem : 
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There are three groups of four routines for the profile solver, fronted solver and 
eigenvalue extraction. The profile solution package consists of: 

SUBROUTINE COMPRO * Sets up the integer array for the profile of global 
stiffness equations to be stored in a profile form. 

SUBROUTINE ASSEM5 - Assembles the element stiffness equations into the 
global equation system stored in a profile form. 

SUBROUTINE SOLUT1 - Controls the iterative solution processes including 
the vector update required for the quasi- and secant-Newton iterations. 

SUBROUTINE DECOMP - Factorizes the global stiffness equations stored in a 
profile form. 

SUBROUTINE SOLVER - Performs the back substitution and generates the 
update vector for the incremental displacement. 

The frontal solution package consists of: 

SUBROUTINE FRONTW - Estimates the front matrix size to be accommodated 
in the core memory. 

SUBROUTINE INTFR - Allocates memory for the work space required for the 
frontal solution. 

SUBROUTINE PRFRNT - Sets up the elimination table for the frontal solution. 

SUBROUTINE FRONTF - Assembles and factorizes die global stiffness equation 
simultaneously. 

SUBROUTINE FRONTB - Performs the back substitution and generates the 
updates for displacement vector. 

SUBROUTINE FRONTR - Controls the iterative solution processes including 
the vector update required for the quasi- and secant-Newton iterations. Calls 
SUBROUTINE FRONTB for the incremental displacement update vector. 

SUBROUTINE VDSKIO - Controls the data stream stored in-core buffer area 
and the actual out-of-core storage devices. 

The eigenvalue analysis package consists of: 

SUBROUTINE EIGENV - Controls the execution of eigenvalue extraction 
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subsystem. 

SUBROUTINE INIMOP - Initializes the array for eigenvectors. 

SUBROUTINE SUBSPC - Performs the subspace iteration and generates a 
specified number of eigenvalues and eigenvectors. 

SUBROUTINE JACOBI - Solves the eigenvalue problem in the subspace by the 
Jacobi iteration. 

There are a number of subprograms used commonly by Algebraic Operation Sub- 
system: 

SUBROUTINE STRUCT - Controls the memory allocation for the global 
algebraic manipulations at the beginning of every analysis. 

SUBROUTINE INITI2 - Allocates memory required for the storage of global 
stiffness matrix and other vectors required in the linear algebraic 
manipulation of finite element equations. 

SUBROUTINE LINESR - Calculates the search distance when the line search 
option is turned on. 

Element Assembly Submodules : 

These are subprograms constructing vectors and matrices appearing in the algo- 
rithmic description of mixed iterative process discussed in die previous subsections. 


SUBROUTINE ASSEM1 - Assembles the displacement stiffness matrix for the 
preconditioning purposes. All the options for kinematic and constitutive 
assumptions are treated on this module. 

SUBROUTINE ASSEM2 - Assembles the coefficient matrix for the transient 
time integration by the Newmark family of algorithm. This routine is 
evolved from SUBROUINE ASSEM1 and handles all the kinematic and 
constitutive options. 

SUBROUTINE ASSEM3 - Assembles the coefficient matrix for the quasi-static 
analysis using the frontal solution subsystem. Large displacement, stress 
stiffening and centrifugal terms are not available in this package. 

SUBROUTINE ASSEM4 - Calculates the nodal strain and recovers the residual 
vector in a mixed form. The subelement solution module is entered from 
this subprogram. 
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Element Loop Structure and Library Routines : 

In the element assembly submodules, element arrays are generated in the loops 
over elements. The protocol to access to the element library is designed and 
implemented which involves a sequence of subroutine calls: 

SUBROUTINE ELVULV - Sets up the current element parameters from the 
element library table.(See Table A1 for variable definition and its values) 

SUBROUTINE CNODEL - Pulls out quantities for the current element from 
the global nodal array and restores them in the element workspace. 
Coordinate transformations necessary for beam and shell elements are 
performed in this subprogram. 

SUBROUTINE DERIV - Sets up the displacement-strain matrix for the current 
element by calling the element library subroutines. Those are: 

SUBROUTINE BPSTRS for plane stress elements, types 3 and 101. 

SUBROUTINE BPSTRN for plane strain elements, types 11 and 102. 

SUBROUTINE BSOLID for three-dimensional solid element, type 7. 

SUBROUTINE BSHELL for three-dimensional shell element, type 75. 

SUBROUTINE BAXSYM for axisymmetric solid-of-revolution elements, 
types 10 and 103. 

SUBROUTINE BTBEAM for linear Timoshenko beam element, type 98. 

SUBROUTINE BASPST for the assumed stress plane stress element, type 

151. 

SUBROUTINE BASPSN for the assumed stress plane strain element, type 

152. 

SUBROUTINE BASSO L for the assumed stress three-dimensional solid 
element, type 154. 

SUBROUTINE UDERIV - A slot for user coded element B matrix routine. 

The following subprograms are used to calculate terms appearing in the finite 
element equations: 
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SUBROUTINE LMPMAS - Calculates nodals weight factor for the strain 
projection. 

SUBROUTINE STIFF - Performs matrix triple products to assemble the element 
stiffness matrix and the element load vector associated with the initial strain 
terms. 

SUBROUTINE STRAIN - Calculates element strain at specified sampling points 
and projects to nodes. 

SUBROUTINE CNSMAS - Assembles the consistent mass matrix for the modal 
and transient analysis. 

SUBROUTINE INITST - Generates initial stress terms for the quasi-static, 
buckling and modal analysis of prestressed structures. 

SUBROUTINE CENMAS - Evaluates the centrifugal mass terms for rotating 
structures at speed. 

SUBROUTINE RESID - Calculates the element residual vector for the global 
element. 

SUBROUTINE SUBFEM - Performs the subelement solution and calculates the 
element residual vector for the subelement mesh. 

SUBROUTINE RELDFG - Calculates relative deformation gradient at the 
element sampling point and projects to node. 

SUBROUTINE RESDYN - Calculates the contribution of mass and damping 
terms in the element residual vector when the transient dynamics option is 
used. 

Material Library : 

A system of subroutines is included in the MHOST code which covers a wide range 
of material models and initial strain assumptions: 

SUBROUTINE SIMPLE - Integrates the stress over the increment assuming the 
elastic-plastic response of the material is represented by total secant modulus 
of elasticity. Also generates the material modulus matrix. 

SUBROUTINE PLASTS - Integrates the stress over the increment and calculates 
the plastic strain by using the radial return algorithm. 


S.Nakazawa 


May, 1991 



MHOST Theoretical Manual 


62 


SUBROUTINE PLASTD - Calculates the consistent elastic-plastic modulus if die 
incremental equivalent plastic strain is positive. 

SUBROUTINE WALKEQ - Integrates the Walker unified creep plasticity 
constitutive equation and also generates a material modulus based on the 
temperature dependent elasticity assumption. 

SUBROUTINE LELAST - Calculates the stress for the constant material modulus 
given as data. 

SUBROUTINE THRSTN - Calculates the thermal strain. 

SUBROUTINE CRPSTN - Calculates the creep strain. 


S.Nakazawa 


May, 1991 


MHOST Theoretical Manual 


63 


REFERENCES 


Abdel Rahman, H. H. (1982), Computational Models for the Nonlinear Analysis of 
Reinforced Concrete Flexural Slab System , Ph.D. Thesis, University College of 
Swansea. 

Annigeri, B. S. (1984), Surface Integral Finite Element Hybrid Method for Localized 
Problems in Continuum Mechanics, Ph.D. Thesis, Department of Mechanical 
Engineering, Massachusetts Institute of Technology. 

Atluri, S. N., P. Tong and H. Murakawa (1983), Recent Studies in Hybrid and Mixed 
Finite Element Methods in Mechanics, Chap. 3 of Hybrid and Mixed Finite Element 
Methods (S.N. Atluri, R.H. Gallagher and O.C. Zienkiewicz, eds.), Wiley, 
Chichester. 

Axelsson, O. (1976), Solution of Linear Systems of Equations: Iterative Methods, In 
'Sparse Matrix Techniques', Lecture Notes in Mathematics, Vol. 572, Springer- 
Verlag, Berlin. 

Axelsson, O. (1978), On Iterative Solution of Large Sparse Systems of Equations 
with Particular Emphasis on Boundary value Problems, TICOM Report 78-4, The 
Texas Institute for Computational Mechanics, University of Texas at Austin, TX. 

Babuska, I., O. C. Zienkiewicz, J. Gago and E. R. de A. Oliveira (eds.) (1986), Accuracy 
Estimates and Adaptive Refinements in Finite Element Computations, Wiley, 
Chichester. 

Babuska, I. (1973), The Finite Element Method with Lagrange Multipliers, Numer. 
Math, 20, 179-192. 

Babuska, I. and A. Miller (1984 A), The Post-Processing Approach in the Finite 
Element Method - Part 1: Calculation of Displacements, Stresses, and Other Higher 
Derivatives of the Displacements, Int. }. num. Meth Eng., 20, 1085-1109. 

Babuska, I. and A. Miller (1984B), The Post-Processing Approach in the Finite 
Element Method - Part 2: The Calculation of Stress Intensity Factors, Int. }. num. 
Meth. Eng., 20, 1111-1129. 


S.Nakazawa 


March, 1991 



MHOST Theoretical Manual 


64 


Bazeley, G. P., Y. K. Cheung, B. M. Irons and O. C. Zienkiewicz (1966), Triangular 
Elements in Plate Bending. Conforming and Non-Conforming Solutions, Proc. 1st 
Conf. Matrix Methods in Structural Mechanics , AFFDL- TR-CC-80, pp 547-576. 

Beer, G. and W. Haas (1982), A Partitioned Frontal Solver for Finite Element 
Analysis, Int. J. num. Meth. Eng., 18, 1623-1654. 

Belytschko, T. (1986), A Review of Recent Developments in Plate and Shell 
Elements, Computational Mechanics - Advances and Trends (ed. A. K. Noor), 
ASME AMD - Vol. 75, pp 217-232. 

Belytschko, T. and W. E. Bachrach (1985), Simple Quadrilaterals with High Coarse 
Mesh Accuracy, Hybrid and Mixed Finite Element Methods (R. L. Spilker and K. W. 
Reed, eds.), ASME, AMD Vol. 73. 

Belytschko,T., C.S.Tsay and W.K.Liu (1981), A Stabilization Matrix for Bilinear 
Mindlin Plate Element, Comp.Meth.Appl.Mech.Eng., 29, 313-327. 

Belytschko, T. and C.-S. Tsay (1983), A Stabilization Procedure for the Quadrilateral 
Plate Element with One-Point Quadrature, lnt. J. num. Meth. Eng., 19, 405-419. 

Benson, D. J. and J. O. Hallquist (1986), A Simple Rigid Body Algorithm for 
Structural Dynamics Programs, lnt. J. Num. Meth. Eng., 22, 723-749. 

Brezzi, F. (1974), On the Existence, Uniqueness and Approximation of Saddle-Point 
Problems Arising from Lagrange Multipliers, RAIRO, 22, 129- 151. 

Brown, P. N. and A. C. Hindmarsh (1986), Matrix-Free Meethods for Stiff System of 
ODE’s, SIAM J. Numer. Anal, 23, 610-638. 

Cantin, G., G. Loubignac and G. Touzot (1978), An Iterative Algorithm to Build 
Continuous Stress and Displacement Solutions, Int. J. num. Meth. Eng., 12, 1493- 
1506. 

Crisfield, M. A. (1975), An Iterative Improvement for Non- Conforming Bending 
Elements, Int. J. num. Meth. Eng., 9, 641-648. 

Crisfield, M. A. (1982), Solution Procedures for Non-Linear Structural Analysis, 
Chap. 1 of Recent Advances in Non-Linear Computational Mechanics (E. Hinton, 
D. R. J. Owen and C. Taylor, eds.), Pineridge Press, Swansea. 

Crisfield, M. A. (1983), An Arc-Length Method Including Line Searches and 
Accelerations, lnt. J. num. Meth. Eng., 19, 1269- 1289. 


S.Nakazawa 


March, 1991 



MHOST Theoretical Manual 


65 


Crisfield,M.(1986), Finite Elements and Solution Procedures for Structural Analysis , 
Vol.l: Linear Analysis, Pineridge Press, Swansea. 

Dias,J.B. and S.Nakazawa (1988), An Approach to Probablistic Finite Element 
Analysis Using a Mixed Iterative Formulation, Paper to be presented at ASME/SES 
Summer Meeting, Unversity of California, Berkeley. 

Fairweather, G. (1978), Finite Element Galerkin methods for Differential Equations, 
Dekker, New York. 

Fine, A.D. (1984), Special Finite Element Displacement Formulations for Modeling 
Singularity Point Behavior, Report PW Engineering /Connecticut Operations. 

Fletcher, C. A. J. (1984), Computastional Galerkin Methods, Springer- Verlag, New 
York. 

Fortin,M. and R.Glowinski(eds.) (1982), Methdes de Lagrangian Augmente, Dunod, 
Paris. 

Gaudreau, G. L., D. J. Benson, J. O. Hallquist, G. J. Kay, R. W. Rosinsky and S. 
J.Sachett (1986), Supercomputers for Engineering Analysis, Computational 
Mechanics - Advances and Trends (ed. A. K, Noor), ASME AMD - Vol. 75, pp 117- 
131. 

Girault, V. and P.-A. Raviart (1979), Finite Element Approximation of the Navier- 
Stokes Equations, Lecture note in Mathematics, 749, Springer-Verlag, Berlin. 

Glowinski, R., E. Y. Rodin and O. C. Zienkiewicz (eds.) (1979), Energy Methods in 
Finite Element Analysis, Wiley, Chichester. 

Haugender, E. and H. A. Mang (1979), On an Improper Modification of a Variational 
Principle for Finite Element Analysis, TAMM, 59, 637-640. 

Hayes, L. J. and P. Devloo (1986), A Vectorized Version of a Spare Matrix- Vector 
Multiply, Int. J. Num. Meth. Eng., 23, 1043-1056. 

Hestens, M. R. (1969), Multiplier and Gradient Methods, J. Optim. Theory Appl., 4, 
303-320. 

Hughes, T. J. R. and T. E. Tezduyar (1981), Finite Elements based upon Mindlin Plate 
Theory with Particular Reference to the Four- Node Bilinear Isoparametric Element, 
J. Appl. Mech. Trans. ASME, 48, 587-596. 


S.Nakazawa 


March, 1991 


MHOST Theoretical Manual 


66 


Hughes, T. J. R. and D. S. Malkus (1983), A General Penalty /Mixed Equivalence 
Theorem for Anisotropic, Incompressible Finite Elements, Chap. 25 of Hybrid and 
Mixed Finite Element Methods (S. N. Atluri, R. H. Gallagher and O. C. Zienkiewicz, 
eds.), Wiley, Chichester. 

Hughes, T. J. R. and W. K. Liu (1981), Nonlinear Finite Element Analysis of Shells; 
Part I: Three Dimensional Shells, Comp. Meth. Appl Mech. Eng., 26, 331-362, 

Hughes, T. J. R. (1983), Analysis of Transient Algorithms with Particular Reference 
to Stability Behavior, Computational Methods for Transient Analysis (eds. T. 
Belytschko and T. J. R. Hughes), Chap. 2, North-Holland, Amsterdam. 

Hughes, T. J. R. and Shakie, F. (1986), Pseudo-Corner Theory: A Simple 
Enhancement of J - Flow Theory for Applications Involving Non-Proportional 
Loading, Eng. Comp., 3, 116-120. 

Hughes, T. J. R., Ferencz, R. M. and Hallquist, J. O. (1986), Large-Scale Vectorized 
Implicit Calculations in Solid Mechanics on a CRAY X-MP/48 Utilizing EBE 
Preconditioned Conjugate Gradients, Computational Mechanics - Advances and 
Trends (ed. A. K. Noor), ASME AMD - Vol. 75, pp 233-279. 

Jameson, A. (1986), Current Status and Future Directions of Computational 
Transonics, Computational Mechanics - Advances and Trends (ed. A. K. Noor), 
ASME AMD - Vol. 75, pp 329-368. 

Kawahara, M. (1978), Steady and Unsteady Finite Element Analysis of 
Incompressible Viscous Fluids, Finite Element in Fluids Vol. 3 ( R. H. Gallagher, O. 
C. Zienkiewicz, J. T. Oden, M. M. Cecci and C. Tayler, eds. ), Wiley, Chichester. 

Kawai, T. and Yoshimura, N. (1969), Analysis of Large Deflection of Plates by Finite 
Element Method, Int. J. Num. Meth Eng., 1 , 123-133. 

Kikuchi, F. and Y. Ando (1972), A New Variational Functional for Finite Element 
Method and Its Application to Plate and Shell Problem, Nucl. Eng. Des., 21, 95-113. 

Lohner, R. (1986), FEM-FCT and Adaptive Refinement Schemes for Strongly 
Unsteady Flows, Numerical Methods for Compressible Flows - Finite Difference, 
Element and Volume Techniques (eds. T. E. Tezduyar and T. J. R. Hughes), ASME 
AMD - Vol. 78, pp 93-114. 

Malkus, D. S. and T. J. R. Hughes (1978), Mixed Finite Element Method - Reduced 
and Selective Integration Technique: A Unification of Concepts, Comp. Meth. 
Appl. Mech. Eng., 15, 63-81. 


S.Nakazawa 


March, 1991 



MHOST Theoretical Manual 


67 


Melosh, R. J. (1963), Basis of Derivation of Matrices for Direct Stiffness Method, 
A1AA /., 1, 1631-1637. 

Muller, A. and T. J. R. Hughes (1984), Augmented Weak Forms and Element-by- 
Element Preconditioners: Efficient Iterative Strategies for Structural Finite 
Elements - A Preliminary Study, Paper Presented at Symposium on Advances and 
Trends in Structures and Dynamics, Arilington, Virginia. 

Muller, A. (1985) Element-By-Element Iterative Procedures in Structural Finite 
Element Analysis, Ph.D. Thesis, Applied Mechanics Division, Stanford University. 

Nagtegaal, J. C. and J. G. Slater (1981), A Simple Non-Compatible Thin Shell 
Element Based on Discrete Kirchhoff Theory, Nonlinear Finite Element Analysis of 
Plates and Shells (eds. T. J. R. Hughes, A. Pifko and A. Jay), ASME AMD - Vol. 48, 
pp 167-192. 

Nagtegaal, J. C., Nakazawa, S. and Tateishi, M. (1986), On the Construction of 
Optimal Mindlin Type Shell Elements, Finite Element Methods for Plate and Shell 
Structures, Vol. 1: Element Technology (ed. T. J. R. Hughes and E. Hinton), 
Pineridge Press, Swansea, pp 348-364, 

Nakazawa, S., J. -P. Vilotte and O. C. Zienkiewicz (1984), Variational Formulations 
and Their Approximations for Incompressible Problems, Paper presented at 5th 
Int'l. Conf. Finite Element in Flow Problems, Austin, Texas. 

Nakazawa, S. (1984A), Reducible and Irreducible Finite Element Forms for Linearly 
Constrained Problems in Mechanics, Proceeding of Fifth Int l. Conf. on Pres. Ves. 
Tech, Vol. Ill, ASME G00242. 

Nakazawa, S. (1984B), Mixed Finite Elements and Iterative Solution Procedures, 
Innovative Methods in Nonlinear Problems (W. K. Liu, T. Belytschko and K. C. 
Park, eds.), Pineridge Press, Swansea, Chap. 21. 

Nakazawa, S., Y. Owa and Q, C. Zienkiewicz (1984), On the Nodal Stress Recovery 
in the Displacement Finite Element Methods, Internal Report, MARC Analysis 
Research Corporation (submitted for publication). 

Nakaza w a,S . ( 1 987) , 3-D Inelastic Analysis Method for Hot Section Components, 
Fourth Annual Status Report: Volume 1: Special Finite Element Models, Pratt & 
Whitney PWA-5940-62. 

Nakazawa,S.(1986), 3-D Inelastic Analysis Method for Hot Section Components, 
Third Annual Status Report: Volume 1: Special Finite Element Models, NASA CR- 
179694, Pratt & Whitney PWA-5940-46. 


S.Nakazawa 


March, 1991 



MHOST Theoretical Manual 


68 


Nakazawa,S.,J.C.Nagtegaal and O.C.Zienkiewicz (1985), Iterative Methods for Mixed 
Finite Element Equations,(R.L.Spilker and K.W.Reed,eds.) Hybrid and Mixed Finite 
Element Methods, ASME AMD Vol.73, pp57-67. 

Nakazawa, S. and Nagtegaal, J. C. (1986), An Efficient Numerical Procedure for 
Nonlinear Analysis of Turbine Engine Hot Section Components, Presented at 3rd 
Symp. on Nonlinear Constitutive Relations for High Temperature Applications, 
Akron, OH. 

Nakazawa, S. and Spiegel, M. S. (1986), Mixed Finite Elements and Iterative 
Solution Processes, Presented at 1st WCCM, Austin, Texas. 

Needleman, A. (1986), Finite Element Analysis of Failure Modes in Ductile Solids, 
Proc. 1st WCCM, Vol. 1, Texas Institute for Computational Mechanics. 

Nour-Omid, B. and M. Ortiz (1986), Concurrent Solution Methods for Finite 
Element Equations, Proc. 1st WCCM, Vol. 1, Texas Institute for Computational 
Mechanics. 

Oden, J. T. and H. J. Brauchli (1971), On the Calculation of Consistent Stress 
Distributions in Finite Element Approximations, lnt. J. num. Meth. Eng., 3, 371-325. 

Oden, J. T. and J. N. Reddy (1973), Note on an Approximate Method for Computing 
Consistent Conjugate Stresses in Elastic Finite Elements, lnt. J. num. Meth. Eng., 6, 
55-61. 


Owa, Y. (1982), On Stress Recovery in Displacement Finite Element Method, M.Sc. 
Thesis, University College of Swansea, Swansea. 

Ortiz, M. and Y. Leroy (1986), A Finite Element Method for Localized Failure 
Analysis, Proc. 1st WCCM, Vol. 1, Texas Institute for Computational Mechanics. 

Ortiz, M. and Simo, J. C. (1986), An Analysis of a New Class of Integration 
Algorithms for Elastoplastic Constitutive Relations, lnt. }. Num. Meth. Eng., 23, 
355-366. 

Peraire, J., K. Morgan and O. C. Zienkiewicz (1986), Convection Dominated 
Problems, Numerical Methods for Compressible Flows - Finite Difference, Element 
and Volume Techniques (eds. T. E. Tezduyar and T. J. R. Hughes), ASME AMD - 
Vol. 78, pp 129-147. 

Pian, T. H. H. and Z. Tian (1985), Hybrid Solid Element with Traction-Free 
Cylindrical Surface, Hybrid and Mixed Finite Element Methods (R. L. Spilker and K. 


S.Nakazawa 


March, 1991 


MHOST Theoretical Manual 


69 


W. Reed, eds.), ASME, AMD Vol. 73. 

Pian, T. H. H. (ed) (1984), Special Memorial Issue Dedicated to Prof. Kyuichiro 
Washizu, Comp. Struct., 19 (1/2). 

Pian,T.H.H., K.Sumihara and D.Kang(1983), New Variational Formulations for 
Hybrid Stress Elements, Nonlinear Structural Analysis (C.C.Chamis,ed.), NASA 
Conference Publication 2297. 

Powell, M. J. D. (1969), A Method for Nonlinear Constraints in Optimization 
Problems (R. Fletcher, ed.). Optimization, Academic Press, London, pp 283 - 298. 

Simo, J.C., R.L. Taylor and K.S. Pister (1984), Variational and Projection Methods 
for the Volume Constraint in Finite Deformation of Elastoplasticity, to appear in 
Comp. Meth. Appl. Mech. Eng. 

Simo, J.C. and T.J.R. Hughes (1985), On the Variational Foundations of Assumed 
Strain Methods, a manuscript submitted for publication. 

Simo, J. C. (1986), A Framework for Finite Strain Elastoplasticity Based on 
Maximum Plastic Dissipation and the Multiplicative Decomposition. Part I: 
Continuum Formulation; Part 11: Computational Aspects, Manuscript, Applied 
Mechanics Division, Stanford University. 

Stlarski, H. and T. Belytschko (1985), Limitation Principles for Mixed Finite 
Elements Based on the Hu - Washizu Variational formulation, Hybrid and Mixed 
Finite Element Methods (R. L. Spilker and K. W. Reed, eds.), ASME, AMD Vol. 73. 

Strang, G. and G. J. Fix (1973), An Analysis of the Finite Element Method, Prentice- 
Hall, Englewood Cliffs, N.J. 

Stummel, F. (1980), The Limitations of Patch Test, Int. J. Num. Meth. Eng., 15, 177- 
188. 

Le Tallec, P. (1986), Mixed Numerical Methods in Viscoplasticity, Proc. 1st WCCM, 
Vol. 1, Texas Institute for Computational Mechanics. 

Taylor, R.L. and O.C. Zienkiewicz (1982), Mixed Finite Element Solution of Fluid 
Flow Problems, Chap. 1 of Finite Elements in Fluids, Vol. 4 (R.M. Gallagher, D.H. 
Norrie, J.T. Oden and O.C. Zienkiewicz, eds.), Wiley, Chichester. 

Taylor, R. L. (1985), Solution of Linear Equations by a Profile Solver, Eng. Comp., 2, 
344-350. 


S.Nakazawa 


March, 1991 



MHOST Theoretical Manual 


70 


Taylor, R. L., Simo, J. C., Zienkiewicz, O. C., and Chen, A. C. H. (1986), The Patch Test 
- A Condition for Assessing FEM Convergence, lnt. J. Num. Meth. Eng., 22, 39-62. 

Temam, R. (19 77), Navier-Stokes Equations, North Holland, Amsterdam. 

Todd, E. S., B. N. Cassenti, S. Nakazawa and P. K. Banerjee (1984), 3-D Inelastic 
Analysis Methods for Hot Section Components (Base Program) First Annual Status 
Report, PWA-5940-19, United Technologies, Pratt and Whitney Aircraft, Hartford, 
CT. 

Utku, S., Salama, M. and Melosh, R. J. (1986), Concurrent Cholesky Factorization of 
Positive Definite Banded Hermitian Matrices, lnt. J. Num. Meth. Eng., 23, 2137- 
2152. 

Utku, S., Melosh, R. J., Salama, M. and Chang, H. T. (1986), Problem Decomposition 
for Parallel Processing, Proc. 1st WCCM, Vol. 1, Texas Institute for Computational 
Mechanics. 

Fraejis de Veubeck, B. (1965), Displacement and Equilibrium Models in the Finite 
Element Method, Stress Analysis ( O. C. Zienkiewicz and G. S. Holister, eds. ), 
Wiley, Chichester. 

Washizu, K. (1955), On the Variational Principles of Elasticity and Plasticity, 
Technical Report 25-18, Aeroelastic and Struct. Res. Lab., MIT. 

Washizu,K.(1974), Variational Methods in Elasticity and Plasticity, Pergamon Press, 
Oxford. 

Way, S. (1938), Uniformly Loaded, Clamped, Rectangular Plates with Large 
Deformation, Proc. 5th Int'l. Cong. Appl. Mech., Cambridge, MA. 

Wiliam, K., Sobh, N. and Sture, S. (1986), Bifurcation Studies of Elastic- Plastic 
Solids, Proc. 1st WCCM, Vol. 1, Texas Institute for Computational Mechanics. 

Wilson,R.B.,M.J.Bak, S.Nakazawa and P.K.Banerjee (1984), 3-D Inelastic Analysis 
Methods for Hot Section Components (Base Program), First Annial Status Report, 
NASA CR- 174700. 

Wilson,R.B.,M.J.Bak, S.Nakazawa and P.K.Banerjee (1986), 3-D Inelastic Analysis 
Methods for Hot Section Components (Base Program), Second Annial Status 
Report, NASA CR-175060. 

Xue, W-M. and S. N. Atluri (1985), Existence and Stability, and Discrete BB and Rank 
Conditions, for General Mixed-Hybrid finite Elements in Elasticity, Hybrid and 


S.Nakazawa 


March, 1991 



MHOST Theoretical Manual 


71 


Mixed Finite Element Methods (R. L. Spilker and K. W. Reed, eds.), ASME, AMD 
Vol. 73. 

Zienkiewicz, O. C. and S. Valliappan (1969), Analysis od Real Structures for Creep, 
Plasticity and Other Complex Constitutive Laws, Proc. Conf. on Civil Engineering 
Materials, Southampton. 

Zienkiewicz, O. C. (1977), The Finite Element Method, 3rd ed., McGraw Hill, New 
York. 

Zienkiewicz, O. C. and S. Nakazawa (1982), The Penalty Function Method and Its 
applications to the Numerical Solution of Boundary Value Problems, AMD Vol. 51, 
ASME, 157-179. 

Zienkiewicz, O. C. and S. Nakazawa (1984), On Variational Formulation and Its 
Modifications for Numerical Solution, Int. J. Solid Struct, (to appear). 

Zienkiewicz, O. C., Li Xi-Kui and S. Nakazawa (1984), Iterative Solution of Mixed 
Problems and the Stress Recovery procedures, Report No. C/R/476/84, Institute for 
Numerical Methods in Engineering, University College of Swansea. 

Zienkiewicz, O. C., J. -P. Vilotte, S. Nakazawa and S. Toyoshima (1984), Iterative 
Method for Constrained and Mixed Approximation. An Inexpensive Improvement 
of F.E.M. Performance, Institute for Numerical Methods in Engineering, Report 
No. C/R/489/84, University College of Swansea (to appear in Comp. Meth. Appl. 
Mech. Eng.) 

Zienkiewicz, O. C., S. Qu, R. L. Taylor and S. Nakazawa (1986), The Patch Test for 
Mixed Formulations, Int. J. Num. Meth. Eng., 23, 1873-1883. 

Zienkiewicz, O. C. and J. Z. Zhu (1987), A Simple Error Estimator and Adaptive 
Procedure for Practical Engineering Analysis, Int. J. Num. Meth. Eng., 24, 337-357. 

Zienkiewicz, O.C., X.K.Li and S.Nakazawa (1984), Iterative Solution of Mixed 
Problems and the Stress Recovery Procedures, Comm. Appl.Numer. Anal, 1, 3-10. 

Zienkie wicz,0 .C ., X.K.Li and S.Nakazawa (1986), Dynamic Transient Analysis by 
Mixed, Iterative Method, IntJ.num.Meth.Eng., 23, 1343-1353. 

Zienkiewicz,O.C. and S.Nakazawa (1984), On Variational Formulation and Its 
Modifications for Numerical Solution, Computer 6- Struct., 19, 303-313. 

Zienkiewicz,O.C.,J P. Vilotte, S.Toyoshima and S.Nakazawa (1984), Iterative Methods 
for Constrained and Mixed Approximation. An Inexpensive Improvement of FEM 


S.Nakazawa 


March, 1991 


MHOST Theoretical Manual 


72 


Performance, Comp.Meth.Appl.Mech.Eng., 51, 3-29. 

Zienkiewicz, O. C. (1986), Three-Field Mixed Approximation and the Plate Bending 
Problem, Report No. C/R/566/86, Institute for Numerical Method in Engineering, 
University College of Swansea. 

Zyvoloski, G. (1986), Incomplete Factorization for Finite Element Methods, Int. J. 
Num. Meth. Eng., 23, 1101-1109. 


S.Nakazawa 


March, 1991 



